diff options
Diffstat (limited to 'src/Skeleton_blocker')
27 files changed, 6171 insertions, 0 deletions
diff --git a/src/Skeleton_blocker/concept/SkeletonBlockerDS.h b/src/Skeleton_blocker/concept/SkeletonBlockerDS.h new file mode 100644 index 00000000..1fa9dc06 --- /dev/null +++ b/src/Skeleton_blocker/concept/SkeletonBlockerDS.h @@ -0,0 +1,127 @@ +/* + * SkeletonBlockerDS.h + * + * Created on: Feb 20, 2014 + * Author: David Salinas + * Copyright 2013 INRIA. All rights reserved + */ + +#ifndef GUDHI_SKELETONBLOCKERDS_H_ +#define GUDHI_SKELETONBLOCKERDS_H_ + + +namespace Gudhi { + +namespace skbl { + + + +/** \brief Concept that must be passed to + * the template class Skeleton_blockers_complex + * + */ +struct SkeletonBlockerDS +{ + /** + * @todo faire un default value pour les vertex_handle + */ + + /** + * @brief index that allows to find the vertex in the boost graph + */ + typedef int boost_vertex_handle; + + + /** + * @brief Root_vertex_handle and Vertex_handle are similar to global and local vertex descriptor + * used in <a href="http://www.boost.org/doc/libs/1_38_0/libs/graph/doc/subgraph.html">boost subgraphs</a> + * and allow to localize a vertex of a subcomplex on its parent root complex. + * + * In gross, vertices are stored in a vector + * and the Root_vertex_handle and Vertex_handle store indices of a vertex in this vector. + * + * For the root simplicial complex, the Root_vertex_handle and Vertex_handle of a vertex + * are the same. + * + * + * For a subcomplex L of a simplicial complex K, the local descriptor, ie the Vertex_handle, of a + * vertex v (that belongs to L) is its position in the vector of vertices + * of the subcomplex L whereas its Root_vertex_handle (global descriptor) is the position of v in the vector of the + * vertices of the root simplicial complex K. + */ + struct Root_vertex_handle{ + + boost_vertex_handle vertex; + + friend ostream& operator << (ostream& o, const Root_vertex_handle & v); + }; + + /** + * A Vertex_handle must be Default Constructible, Assignable and Equality Comparable. + */ + struct Vertex_handle{ + boost_vertex_handle vertex; + + friend ostream& operator << (ostream& o, const Vertex_handle & v); + }; + + + /** + * \brief The type of vertices that are stored the boost graph. + * A Vertex must be Default Constructible and Equality Comparable. + * + */ + struct Graph_vertex{ + /** \brief Used to deactivate a vertex for example when contracting an edge. + * It allows in some cases to remove the vertex at low cost. + */ + void deactivate(); + + /** \brief Used to activate a vertex. + */ + void activate(); + + /** \brief Tells if the vertex is active. + */ + bool is_active() const; + + void set_id(Root_vertex_handle i); + Root_vertex_handle get_id() const; + virtual string to_string() const ; + friend ostream& operator << (ostream& o, const Graph_vertex & v); + }; + + + /** + * \brief The type of edges that are stored the boost graph. + * An Edge must be Default Constructible and Equality Comparable. + */ + struct Graph_edge{ + /** + * @brief Allows to modify vertices of the edge. + */ + void setId(Root_vertex_handle a,Root_vertex_handle b); + + /** + * @brief Returns the first vertex of the edge. + */ + Root_vertex_handle first() const ; + + /** + * @brief Returns the second vertex of the edge. + */ + Root_vertex_handle second() const ; + + friend ostream& operator << (ostream& o, const Simple_edge & v); + }; + + +}; + + + +} // namespace skbl +} // namespace GUDHI + + +#endif /* GUDHI_SKELETONBLOCKERDS_H_ */ diff --git a/src/Skeleton_blocker/concept/SkeletonBlockerGeometricDS.h b/src/Skeleton_blocker/concept/SkeletonBlockerGeometricDS.h new file mode 100644 index 00000000..23edaaef --- /dev/null +++ b/src/Skeleton_blocker/concept/SkeletonBlockerGeometricDS.h @@ -0,0 +1,65 @@ +/* + * SkeletonBlockerGeometricDS.h + * + * Created on: Feb 20, 2014 + * Author: David Salinas + * Copyright 2013 INRIA. All rights reserved + */ + +#ifndef GUDHI_SKELETONBLOCKERGEOMETRICDS_H_ +#define GUDHI_SKELETONBLOCKERGEOMETRICDS_H_ + +/** \brief Concept that must be passed to + * the template class Skeleton_blocker_geometric_complex + * + */ +template<typename GeometryTrait> +struct SkeletonBlockerGeometricDS : public SkeletonBlockerDS +{ + + /** + * Geometry information. + */ + typedef GeometryTrait GT ; + + /** + * Type of point (should be the same as GT::Point). + */ + typedef typename GeometryTrait::Point Point; + + /** + * @brief Vertex that stores a point. + */ + class Graph_vertex : public SkeletonBlockerDS::Graph_vertex{ + public: + /** + * @brief Access to the point. + */ + Point& point(){ return point_; } + /** + * @brief Access to the point. + */ + const Point& point() const { return point_; } + }; + + /** + * @brief Edge that allows to access to an index. + * The indices of the edges are used to store heap information + * in the edge contraction algorithm. + */ + class Graph_Edge : public SkeletonBlockerDS::Graph_edge{ + public: + /** + * @brief Access to the index. + */ + int& index(); + /** + * @brief Access to the index. + */ + int index(); + }; +}; + + + +#endif /* GUDHI_SKELETONBLOCKERGEOMETRICDS_H_ */ diff --git a/src/Skeleton_blocker/doc/SkeletonBlockerMainPage.h b/src/Skeleton_blocker/doc/SkeletonBlockerMainPage.h new file mode 100644 index 00000000..5a80bc8d --- /dev/null +++ b/src/Skeleton_blocker/doc/SkeletonBlockerMainPage.h @@ -0,0 +1,156 @@ + +/*! \mainpage Skeleton blockers data-structure + +\author David Salinas + +\section Introduction +The Skeleton-blocker data-structure had been introduced in the two papers +\cite skbl_socg2011 \cite skbl_ijcga2012. +It proposes a light encoding for simplicial complexes by storing only an *implicit* representation of its +simplices. +Intuitively, it just stores the 1-skeleton of a simplicial complex with a graph and the set of its "missing faces" that +is very small in practice (see next section for a formal definition). +This data-structure handles every classical operations used for simplicial complexes such as + as simplex enumeration or simplex removal but operations that are particularly efficient + are operations that do not require simplex enumeration such as edge iteration, link computation or simplex contraction. + + +\todo{image wont include} + +\section Definitions + +We recall briefly classical definitions of simplicial complexes \cite Munkres. +An abstract simplex is a finite non-empty set and its dimension is its number of elements minus 1. +Whenever \f$\tau \subset \sigma\f$ and \f$\tau \neq \emptyset \f$, \f$ \tau \f$ is called a face of +\f$ \sigma\f$ and \f$ \sigma\f$ is called a coface of \f$ \tau \f$ . Furthermore, +when \f$ \tau \neq \sigma\f$ we say that \f$ \tau\f$ is a proper-face of \f$ \sigma\f$. +An abstract simplicial complex is a set of simplices that contains all the faces of its simplices. +The 1-skeleton of a simplicial complex (or its graph) consists of its elements of dimension lower than 2. + + +\image latex "images/ds_representation.eps" "My application" width=10cm + + +To encode, a simplicial complex, one can encodes all its simplices. +In case when this number gets too large, +a lighter and implicit version consists of encoding only its graph plus some elements called missing faces or blockers. +A blocker is a simplex of dimension greater than 1 +that does not belong to the complex but whose all proper faces does. + + +Remark that for a clique complex (i.e. a simplicial complex whose simplices are cliques of its graph), the set of blockers +is empty and the data-structure is then particularly sparse. +One famous example of clique-complex is the Rips complex which is intensively used +in topological data-analysis. +In practice, the set of blockers of a simplicial complex +remains also small when simplifying a Rips complex with edge contractions +but also for most of the simplicial complexes used in topological data-analysis such as Delaunay, Cech or Witness complexes. +For instance, the numbers of blockers is depicted for a random 3 dimensional sphere embedded into \f$R^4\f$ +in figure X. + + +\image latex "images/blockers_curve.eps" "My application" width=10cm + + + + +\section API + +\subsection Overview + +Four classes define a simplicial complex namely : + +\li <Code>Skeleton_blocker_complex</Code> : a simplicial complex with basic operations such as vertex/edge/simplex enumeration and construction +\li <Code>Skeleton_blocker_link_complex</Code> : the link of a simplex in a parent complex. It is represented as a sub complex +of the parent complex +\li <Code>Skeleton_blocker_simplifiable_complex</Code> : a simplicial complex with simplification operations such as edge contraction or simplex collapse +\li <Code>Skeleton_blocker_geometric_complex</Code> : a simplicial complex who has access to geometric points in \f$R^d\f$ + +The two last classes are derived classes from the <Code>Skeleton_blocker_complex</Code> class. The class <Code>Skeleton_blocker_link_complex</Code> inheritates from a template passed parameter +that may be either <Code>Skeleton_blocker_complex</Code> or <Code>Skeleton_blocker_geometric_complex</Code> (a link may store points coordinates or not). + +\todo{include links} + +\subsection Visitor + +The class <Code>Skeleton_blocker_complex</Code> has a visitor that is called when usual operations such as adding an edge or remove a vertex are called. +You may want to use this visitor to compute statistics or to update another data-structure (for instance this visitor is heavily used in the +<Code>Contraction</Code> package). + + + +\section Example + + +\subsection s Iterating through vertices, edges, blockers and simplices + +Iteration through vertices, edges, simplices or blockers is straightforward with c++11 for range loops. +Note that simplex iteration with this implicit data-structure just takes +a few more time compared to iteration via an explicit representation +such as the Simplex Tree. The following example computes the Euler Characteristic +of a simplicial complex. + + \code{.cpp} + typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex; + typedef Complex::Vertex_handle Vertex_handle; + typedef Complex::Simplex_handle Simplex; + + const int n = 15; + + // build a full complex with 10 vertices and 2^n-1 simplices + Complex complex; + for(int i=0;i<n;i++) + complex.add_vertex(); + for(int i=0;i<n;i++) + for(int j=0;j<i;j++) + //note that add_edge adds the edge and all its cofaces + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + + // this is just to illustrate iterators, to count number of vertices + // or edges, complex.num_vertices() and complex.num_edges() are + // more appropriated! + unsigned num_vertices = 0; + for(auto v : complex.vertex_range()){ + ++num_vertices; + } + + unsigned num_edges = 0; + for(auto e : complex.edge_range()) + ++num_edges; + + unsigned euler = 0; + unsigned num_simplices = 0; + // we use a reference to a simplex instead of a copy + // value here because a simplex is a set of integers + // and copying it cost time + for(const Simplex & s : complex.simplex_range()){ + ++num_simplices; + if(s.dimension()%2 == 0) + euler += 1; + else + euler -= 1; + } + std::cout << "Saw "<<num_vertices<<" vertices, "<<num_edges<<" edges and "<<num_simplices<<" simplices"<<std::endl; + std::cout << "The Euler Characteristic is "<<euler<<std::endl; + \endcode + + +\verbatim +./SkeletonBlockerIteration +Saw 15 vertices, 105 edges and 32767 simplices +The Euler Characteristic is 1 + 0.537302s wall, 0.530000s user + 0.000000s system = 0.530000s CPU (98.6%) +\endverbatim + + + +\subsection Acknowledgements +The author wishes to thank Dominique Attali and André Lieutier for +their collaboration to write the two initial papers about this data-structure + and also Dominique for leaving him use a prototype. + + +\copyright GNU General Public License v3. +\verbatim Contact: David Salinas, david.salinas@inria.fr \endverbatim + +*/ diff --git a/src/Skeleton_blocker/example/CMakeLists.txt b/src/Skeleton_blocker/example/CMakeLists.txt new file mode 100644 index 00000000..0aed1e7a --- /dev/null +++ b/src/Skeleton_blocker/example/CMakeLists.txt @@ -0,0 +1,7 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIskbl) + +add_executable ( SkeletonBlockerIteration Skeleton_blocker_iteration.cpp ) + +target_link_libraries( SkeletonBlockerIteration ${Boost_TIMER_LIBRARY}) + diff --git a/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp b/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp new file mode 100644 index 00000000..3f4a9ffc --- /dev/null +++ b/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp @@ -0,0 +1,73 @@ +#include <boost/timer/timer.hpp> +#include <stdio.h> +#include <stdlib.h> +#include <string> +#include <fstream> +#include <sstream> + + +#include "gudhi/Skeleton_blocker.h" +//#include "gudhi/Skeleton_blocker_complex.h" +//#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h" + +using namespace std; +using namespace Gudhi; +using namespace skbl; + +typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex; +typedef Complex::Vertex_handle Vertex_handle; +typedef Complex::Simplex_handle Simplex; + + +Complex build_complete_complex(int n){ + // build a full complex with 10 vertices and 2^n-1 simplices + Complex complex; + for(int i=0;i<n;i++) + complex.add_vertex(); + for(int i=0;i<n;i++) + for(int j=0;j<i;j++) + //note that add_edge, add the edge and all its cofaces + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + return complex; +} + +int main (int argc, char *argv[]){ + boost::timer::auto_cpu_timer t; + + const int n = 15; + + // build a full complex with 10 vertices and 2^n-1 simplices + Complex complex(build_complete_complex(n)); + + // this is just to illustrate iterators, to count number of vertices + // or edges, complex.num_vertices() and complex.num_edges() are + // more appropriated! + unsigned num_vertices = 0; + for(auto v : complex.vertex_range()){ + if(v==v); //do something with v as removing a ennoying warning. + ++num_vertices; + } + + unsigned num_edges = 0; + for(auto e : complex.edge_range()){ + if(e==e); + ++num_edges; + } + + unsigned euler = 0; + unsigned num_simplices = 0; + // we use a reference to a simplex instead of a copy + // value here because a simplex is a set of integers + // and copying it cost time + for(const Simplex & s : complex.simplex_range()){ + ++num_simplices; + if(s.dimension()%2 == 0) + euler += 1; + else + euler -= 1; + } + std::cout << "Saw "<<num_vertices<<" vertices, "<<num_edges<<" edges and "<<num_simplices<<" simplices"<<std::endl; + std::cout << "The Euler Characteristic is "<<euler<<std::endl; + return EXIT_SUCCESS; +} + diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h new file mode 100644 index 00000000..1d215984 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h @@ -0,0 +1,24 @@ +/* + * Skeleton_blocker.h + * + * Created on: Jan, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKER_H +#define GUDHI_SKELETON_BLOCKER_H + +#include "gudhi/Skeleton_blocker_complex.h" +#include "gudhi/Skeleton_blocker_geometric_complex.h" +#include "gudhi/Skeleton_blocker_simplifiable_complex.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h" + +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h" + + +#include "gudhi/Utils.h" //xxx + + +#endif + diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h new file mode 100644 index 00000000..4475afcc --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h @@ -0,0 +1,113 @@ +/* + * ComplexVisitor.h + * + * Created on: Dec 11, 2013 + * Author: dsalinas + */ + +#ifndef GUDHI_COMPLEXVISITOR_H_ +#define GUDHI_COMPLEXVISITOR_H_ + + +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h" + +namespace Gudhi{ + +namespace skbl { +// todo rajouter les const + +/** + *@class Skeleton_blocker_complex_visitor + *@brief Interface for a visitor of a simplicial complex. + */ +template <typename Vertex_handle> +class Skeleton_blocker_complex_visitor { +public: + virtual ~Skeleton_blocker_complex_visitor(){}; + + virtual void on_add_vertex(Vertex_handle) = 0; + virtual void on_remove_vertex(Vertex_handle) = 0; + + virtual void on_add_edge(Vertex_handle a,Vertex_handle b) = 0; + virtual void on_remove_edge(Vertex_handle a,Vertex_handle b) = 0; + + /** + * @brief Called when an edge changes, during the contraction of + * an edge + */ + virtual void on_changed_edge(Vertex_handle a,Vertex_handle b) = 0; + + /** + * @brief Called when performing an edge contraction when + * an edge bx is replaced by an edge ax (not already present). + * Precisely, this methods is called this way in contract_edge : + * add_edge(a,x) + * on_swaped_edge(a,b,x) + * remove_edge(b,x) + */ + virtual void on_swaped_edge(Vertex_handle a,Vertex_handle b,Vertex_handle x)=0; + virtual void on_add_blocker(const Skeleton_blocker_simplex<Vertex_handle>&) = 0; + virtual void on_delete_blocker(const Skeleton_blocker_simplex<Vertex_handle>*) = 0; +}; + + +/** + *@class Dummy_complex_visitor + *@brief A dummy visitor of a simplicial complex that does nothing + * + */ +template <typename Vertex_handle> +class Dummy_complex_visitor : public Skeleton_blocker_complex_visitor<Vertex_handle>{ +public: + void on_add_vertex(Vertex_handle) {} + void on_remove_vertex(Vertex_handle){} + void on_add_edge(Vertex_handle a,Vertex_handle b){} + void on_remove_edge(Vertex_handle a,Vertex_handle b){} + void on_changed_edge(Vertex_handle a,Vertex_handle b){} + void on_swaped_edge(Vertex_handle a,Vertex_handle b,Vertex_handle x){} + void on_add_blocker(const Skeleton_blocker_simplex<Vertex_handle>&){} + void on_delete_blocker(const Skeleton_blocker_simplex<Vertex_handle>*){} + +}; + + + +/** + *@class Print_complex_visitor + *@brief A visitor that prints the visit information. + * + */ +template <typename Vertex_handle> +class Print_complex_visitor : public Skeleton_blocker_complex_visitor<Vertex_handle>{ +public: + void on_add_vertex(Vertex_handle v) { + std::cerr << "on_add_vertex:"<<v<<std::endl; + } + void on_remove_vertex(Vertex_handle v){ + std::cerr << "on_remove_vertex:"<<v<<std::endl; + } + void on_add_edge(Vertex_handle a,Vertex_handle b){ + std::cerr << "on_add_edge:"<<a<<","<<b<<std::endl; + } + void on_remove_edge(Vertex_handle a,Vertex_handle b){ + std::cerr << "on_remove_edge:"<<a<<","<<b<<std::endl; + } + void on_changed_edge(Vertex_handle a,Vertex_handle b){ + std::cerr << "on_changed_edge:"<<a<<","<<b<<std::endl; + } + void on_swaped_edge(Vertex_handle a,Vertex_handle b,Vertex_handle x){ + std::cerr << "on_swaped_edge:"<<a<<","<<b<<","<<x<<std::endl; + } + void on_add_blocker(const Skeleton_blocker_simplex<Vertex_handle>& b){ + std::cerr << "on_add_blocker:"<<b<<std::endl; + } + void on_delete_blocker(const Skeleton_blocker_simplex<Vertex_handle>* b){ + std::cerr << "on_delete_blocker:"<<b<<std::endl; + } +}; + +} + +} // namespace GUDHI + +#endif /* GUDHI_COMPLEXVISITOR_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h new file mode 100644 index 00000000..864ee6ac --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h @@ -0,0 +1,68 @@ +/* + * Skeleton_blocker_link_superior.h + * + * Created on: Feb 19, 2014 + * Author: David Salinas + * Copyright 2013 INRIA. All rights reserved + */ + +#ifndef GUDHI_SKELETON_BLOCKER_LINK_SUPERIOR_H_ +#define GUDHI_SKELETON_BLOCKER_LINK_SUPERIOR_H_ + +#include "gudhi/Skeleton_blocker_link_complex.h" + +namespace Gudhi{ + +namespace skbl { + +template<class ComplexType> class Skeleton_blocker_sub_complex; + + +/** + * \brief Class representing the link of a simplicial complex encoded by a skeleton/blockers pair. + * It computes only vertices greater than the simplex used to build the link. + */ +template<typename ComplexType> +class Skeleton_blocker_link_superior : public Skeleton_blocker_link_complex<ComplexType> +{ + typedef typename ComplexType::Edge_handle Edge_handle; + + typedef typename ComplexType::boost_vertex_handle boost_vertex_handle; + +public: + + typedef typename ComplexType::Vertex_handle Vertex_handle; + typedef typename ComplexType::Root_vertex_handle Root_vertex_handle; + typedef typename ComplexType::Simplex_handle Simplex_handle; + typedef typename ComplexType::Root_simplex_handle Root_simplex_handle; + typedef typename ComplexType::BlockerMap BlockerMap; + typedef typename ComplexType::BlockerPair BlockerPair; + typedef typename ComplexType::BlockerMapIterator BlockerMapIterator; + typedef typename ComplexType::BlockerMapConstIterator BlockerMapConstIterator; + typedef typename ComplexType::Simplex_handle::Simplex_vertex_const_iterator AddressSimplexConstIterator; + typedef typename ComplexType::Root_simplex_handle::Simplex_vertex_const_iterator IdSimplexConstIterator; + + + Skeleton_blocker_link_superior() + :Skeleton_blocker_link_complex<ComplexType>(true) + { + } + + Skeleton_blocker_link_superior(const ComplexType & parent_complex, Simplex_handle& alpha_parent_adress) + :Skeleton_blocker_link_complex<ComplexType>(parent_complex,alpha_parent_adress,true) + { + } + + Skeleton_blocker_link_superior(const ComplexType & parent_complex, Vertex_handle a_parent_adress) + :Skeleton_blocker_link_complex<ComplexType>(parent_complex,a_parent_adress,true) + { + } + +}; + +} + +} // namespace GUDHI + + +#endif /* SKELETON_BLOCKER_LINK_SUPERIOR_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h new file mode 100644 index 00000000..369635e5 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h @@ -0,0 +1,108 @@ +/* + * Skeleton_blocker_off_io.h + * Created on: Nov 28, 2014 + * 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): David Salinas + * + * Copyright (C) 2014 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 SKELETON_BLOCKER_OFF_IO_H_ +#define SKELETON_BLOCKER_OFF_IO_H_ + +#include "gudhi/Off_reader.h" + +namespace Gudhi { + +namespace skbl { + +template<typename Complex> +class Skeleton_blocker_off_visitor_reader{ + Complex& complex_; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Point Point; + + const bool load_only_points_; + +public: + Skeleton_blocker_off_visitor_reader(Complex& complex,bool load_only_points = false): + complex_(complex), + load_only_points_(load_only_points){} + + + void init(int dim,int num_vertices,int num_faces,int num_edges){ + //todo do an assert to check that this number are correctly read + //todo reserve size for vector points + } + + + void point(const std::vector<double>& point){ + complex_.add_vertex(point); + } + + void maximal_face(const std::vector<int>& face){ + if (!load_only_points_){ + for(size_t i = 0 ; i < face.size();++i) + for(size_t j = i+1 ; j < face.size();++j){ + complex_.add_edge(Vertex_handle(face[i]),Vertex_handle(face[j])); + } + } + } + + void done(){ + } +}; + +template<typename Complex> +class Skeleton_blocker_off_reader{ +public: + /** + * name_file : file to read + * read_complex : complex that will receive the file content + * read_only_points : specify true if only the points must be read + */ + Skeleton_blocker_off_reader(const std::string & name_file,Complex& read_complex,bool read_only_points = false):valid_(false){ + std::ifstream stream(name_file); + if(stream.is_open()){ + Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex,read_only_points); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } + } + + /** + * return true iff reading did not meet problems. + */ + bool is_valid() const{ + return valid_; + } + +private: + bool valid_; +}; + +} // namespace skbl + + +} // namespace Gudhi + + +#endif /* SKELETON_BLOCKER_OFF_IO_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h new file mode 100644 index 00000000..79a8e7aa --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h @@ -0,0 +1,66 @@ +/* + * Skeleton_blocker_simple_geometric_traits.h + * + * Created on: Feb 11, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_SIMPLE_GEOMETRIC_TRAITS_H_ +#define GUDHI_SKELETON_BLOCKERS_SIMPLE_GEOMETRIC_TRAITS_H_ + +#include <string> +#include <sstream> +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h" + +namespace Gudhi{ + +namespace skbl{ + +/** + * @extends SkeletonBlockerGeometricDS + */ +template<typename GeometryTrait> +struct Skeleton_blocker_simple_geometric_traits : public skbl::Skeleton_blocker_simple_traits { +public: + + typedef GeometryTrait GT; + typedef typename GT::Point Point; + typedef typename Skeleton_blocker_simple_traits::Root_vertex_handle Root_vertex_handle; + typedef typename Skeleton_blocker_simple_traits::Graph_vertex Simple_vertex; + + /** + * @brief Vertex with a point attached. + */ + class Simple_geometric_vertex : public Simple_vertex{ + template<class ComplexGeometricTraits> friend class Skeleton_blocker_geometric_complex; + private: + Point point_; + Point& point(){ return point_; } + const Point& point() const { return point_; } + }; + + + class Simple_geometric_edge : public Skeleton_blocker_simple_traits::Graph_edge{ + int index_; + public: + Simple_geometric_edge():Skeleton_blocker_simple_traits::Graph_edge(),index_(-1){} + /** + * @brief Allows to modify the index of the edge. + * The indices of the edge are used to store heap information + * in the edge contraction algorithm. + */ + int& index(){return index_;} + int index() const {return index_;} + }; + + + typedef Simple_geometric_vertex Graph_vertex; + typedef Skeleton_blocker_simple_traits::Graph_edge Graph_edge; +}; + + +} + +} // namespace GUDHI + +#endif /* GUDHI_SKELETON_BLOCKERS_SIMPLE_GEOMETRIC_TRAITS_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h new file mode 100644 index 00000000..3975bb46 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h @@ -0,0 +1,142 @@ +/* + * Skeleton_blocker_simple_traits.h + * + * Created on: Feb 11, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_SIMPLE_TRAITS_H_ +#define GUDHI_SKELETON_BLOCKERS_SIMPLE_TRAITS_H_ + +#include <string> +#include <sstream> +#include "Skeleton_blocker_simplex.h" + +namespace Gudhi{ + +namespace skbl { + +/** + * @extends SkeletonBlockerDS + */ +struct Skeleton_blocker_simple_traits{ + /** + * @brief global and local handle similar to <a href="http://www.boost.org/doc/libs/1_38_0/libs/graph/doc/subgraph.html">boost subgraphs</a>. + * In gross, vertices are stored in a vector. + * For the root simplicial complex, the local and global descriptors are the same. + * For a subcomplex L and one of its vertices 'v', the local descriptor of 'v' is its position in + * the vertex vector of the subcomplex L whereas its global descriptor is the position of 'v' + * in the vertex vector of the root simplicial complex. + */ + struct Root_vertex_handle{ + typedef int boost_vertex_handle; + explicit Root_vertex_handle(boost_vertex_handle val = -1 ):vertex(val){} + boost_vertex_handle vertex; + + bool operator!=( const Root_vertex_handle& other) const{ + return ! (this->vertex == other.vertex); + } + + bool operator==( const Root_vertex_handle& other) const{ + return this->vertex == other.vertex; + } + + bool operator<( const Root_vertex_handle& other) const{ + return this->vertex < other.vertex; + } + + friend std::ostream& operator << (std::ostream& o, const Root_vertex_handle & v) + { + o << v.vertex; + return o; + } + }; + + struct Vertex_handle{ + typedef int boost_vertex_handle; + Vertex_handle(boost_vertex_handle val=-1):vertex(val){} + + boost_vertex_handle vertex; + + bool operator==( const Vertex_handle& other) const{ + return this->vertex == other.vertex; + } + + bool operator!=( const Vertex_handle& other) const{ + return this->vertex != other.vertex; + } + + bool operator<(const Vertex_handle& other) const{ + return this->vertex < other.vertex; + } + + friend std::ostream& operator << (std::ostream& o, const Vertex_handle & v) + { + o << v.vertex; + return o; + } + }; + + + + class Graph_vertex { + bool is_active_; + Root_vertex_handle id_; + public: + virtual ~Graph_vertex(){} + + void activate(){is_active_=true;} + void deactivate(){is_active_=false;} + bool is_active() const{return is_active_;} + void set_id(Root_vertex_handle i){id_=i;} + Root_vertex_handle get_id() const{return id_;} + + virtual std::string to_string() const { + std::ostringstream res; + res<< id_; + return res.str(); + } + + friend std::ostream& operator << (std::ostream& o, const Graph_vertex & v){ + o << v.to_string(); + return o; + } + }; + + class Graph_edge { + Root_vertex_handle a_; + Root_vertex_handle b_; + int index_; + public: + + Graph_edge():a_(-1),b_(-1),index_(-1){} + + int& index(){return index_;} + int index() const {return index_;} + + void setId(Root_vertex_handle a,Root_vertex_handle b){ + a_ = a; + b_ = b; + } + + Root_vertex_handle first() const { + return a_; + } + + Root_vertex_handle second() const { + return b_; + } + + friend std::ostream& operator << (std::ostream& o, const Graph_edge & v){ + o << "("<<v.a_<<","<<v.b_<<" - id = "<<v.index(); + return o; + } + }; + +}; + +} + +} // namespace GUDHI + +#endif /* GUDHI_SKELETON_BLOCKERS_SIMPLE_TRAITS_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h new file mode 100644 index 00000000..5593a5a8 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h @@ -0,0 +1,375 @@ +#ifndef GUDHI_SKELETON_BLOCKER_SIMPLEX_H +#define GUDHI_SKELETON_BLOCKER_SIMPLEX_H + +#include<cassert> +#include<iostream> +#include<set> +#include<vector> +#include <initializer_list> + +namespace Gudhi{ + +namespace skbl { + +/** + *@brief Abstract simplex used in Skeleton blockers data-structure. + * + * An abstract simplex is represented as an ordered set of T elements, + * each element representing a vertex. + * + * The element representing a vertex can be SkeletonBlockerDS::Vertex_handle or SkeletonBlockerDS::Root_vertex_handle. + * + * + */ +template <typename T> + +class Skeleton_blocker_simplex { + +private : + std::set<T> simplex_set; + +public: + typedef typename T::boost_vertex_handle boost_vertex_handle; + + typedef T Vertex_handle; + + + typedef typename std::set<T>::const_iterator Simplex_vertex_const_iterator; + typedef typename std::set<T>::iterator Simplex_vertex_iterator; + + /** @name Constructors / Destructors / Initialization + */ + //@{ + + +// Skeleton_blocker_simplex():simplex_set() {} + + /** + * Clear the simplex + */ + void clear() { + simplex_set.clear(); + } + + + Skeleton_blocker_simplex(std::initializer_list<T>& list) { + for_each(list.begin(),list.end(),add_vertex); + } + + + template<typename... Args> + Skeleton_blocker_simplex(Args... args){ + add_vertices(args...); + } + + + template<typename... Args> + void add_vertices(T v,Args... args){ + add_vertex(v); + add_vertices(args...); + } + + void add_vertices(T v){ + add_vertex(v); + } + + void add_vertices(){ + } + + /** + * Initialize a simplex with a string such as {0,1,2} + */ + Skeleton_blocker_simplex(std::string token){ + clear(); + if ((token[0] == '{') && (token[token.size()-1] == '}' ) ) + { + token.erase(0,1); + token.erase(token.size()-1,1); + while(token.size()!=0 ){ + int coma_position=token.find_first_of(','); + //cout << "coma_position:"<<coma_position<<endl; + std::string n = token.substr(0,coma_position); + //cout << "token:"<<token<<endl; + token.erase(0,n.size()+1); + add_vertex((T)(atoi(n.c_str()))); + } + } + } + + //@} + + /** @name Simplex manipulation + */ + //@{ + + /** + * Add the vertex v to the simplex: + * Note that adding two times the same vertex is + * the same that adding it once. + * \f$ (*this) \leftarrow (*this) \cup \{ v \} \f$ + */ + void add_vertex(T v) + { + simplex_set.insert(v); + } + + /** + * Remove the vertex v from the simplex: + * \f$ (*this) \leftarrow (*this) \setminus \{ v \} \f$ + */ + void remove_vertex(T v) + { + simplex_set.erase(v); + } + + /** + * Intersects the simplex with the simplex a: + * \f$ (*this) \leftarrow (*this) \cap a \f$ + */ + void intersection(const Skeleton_blocker_simplex & a){ + std::vector<T> v; + v.reserve(std::min(simplex_set.size(), a.simplex_set.size())); + + set_intersection(simplex_set.begin(),simplex_set.end(), + a.simplex_set.begin(),a.simplex_set.end(), + std::back_inserter(v)); + clear(); + for (auto i:v) + simplex_set.insert(i); + } + + /** + * Substracts a from the simplex: + * \f$ (*this) \leftarrow (*this) \setminus a \f$ + */ + void difference(const Skeleton_blocker_simplex & a){ + std::vector<T> v; + v.reserve(simplex_set.size()); + + set_difference(simplex_set.begin(),simplex_set.end(), + a.simplex_set.begin(),a.simplex_set.end(), + std::back_inserter(v)); + + clear(); + for (auto i:v) + simplex_set.insert(i); + } + + /** + * Add vertices of a to the simplex: + * \f$ (*this) \leftarrow (*this) \cup a \f$ + */ + void union_vertices(const Skeleton_blocker_simplex & a){ + std::vector<T> v; + v.reserve(simplex_set.size() + a.simplex_set.size()); + + set_union(simplex_set.begin(),simplex_set.end(), + a.simplex_set.begin(),a.simplex_set.end(), + std::back_inserter(v)); + clear(); + simplex_set.insert(v.begin(),v.end()); + } + + typename std::set<T>::const_iterator begin() const{ + return simplex_set.cbegin(); + } + + typename std::set<T>::const_iterator end() const{ + return simplex_set.cend(); + } + + typename std::set<T>::iterator begin(){ + return simplex_set.begin(); + } + + typename std::set<T>::iterator end(){ + return simplex_set.end(); + } + + + + //@} + + /** @name Queries + */ + //@{ + + /** + * Returns the dimension of the simplex. + */ + int dimension() const{ + return (simplex_set.size() - 1); + } + + bool empty() const{ + return simplex_set.empty(); + } + + /** + * Returns the first vertex of the (oriented) simplex. + * + * Be careful : assumes the simplex is non-empty. + */ + T first_vertex() const{ + assert(!empty()); + return *(begin()); + } + + /** + * Returns the last vertex of the (oriented) simplex. + * + * Be careful : assumes the simplex is non-empty. + */ + T last_vertex() const + { + assert(!empty()); + return *(simplex_set.rbegin()); + } + /** + * @return true iff the simplex contains the simplex a, i.e. iff \f$ a \subset (*this) \f$. + */ + bool contains(const Skeleton_blocker_simplex & a) const{ + return includes(simplex_set.cbegin(),simplex_set.cend(),a.simplex_set.cbegin(),a.simplex_set.cend()); + } + + /** + * @return true iff the simplex contains the difference \f$ a \setminus b \f$. + */ + bool contains_difference(const Skeleton_blocker_simplex& a, const Skeleton_blocker_simplex& b) const{ + auto first1 = begin(); + auto last1 = end(); + + auto first2 = a.begin(); + auto last2 = a.end(); + + while (first2!=last2) { + // we ignore vertices of b + if(b.contains(*first2)){ + ++first2; + } + else{ + if ( (first1==last1) || (*first2<*first1) ) return false; + if (!(*first1<*first2)) ++first2; + ++first1; + } + } + return true; + } + + /** + * @return true iff the simplex contains the difference \f$ a \setminus \{ x \} \f$. + */ + bool contains_difference(const Skeleton_blocker_simplex& a, T x) const{ + auto first1 = begin(); + auto last1 = end(); + + auto first2 = a.begin(); + auto last2 = a.end(); + + while (first2!=last2) { + // we ignore vertices of b + if(x == *first2){ + ++first2; + } + else{ + if ( (first1==last1) || (*first2<*first1) ) return false; + if (!(*first1<*first2)) ++first2; + ++first1; + } + } + return true; + } + + /** + * @return true iff the simplex contains the difference \f$ a \setminus \{ x \} \f$. + */ + bool contains_difference(const Skeleton_blocker_simplex& a, T x, T y) const{ + auto first1 = begin(); + auto last1 = end(); + + auto first2 = a.begin(); + auto last2 = a.end(); + + while (first2!=last2) { + // we ignore vertices of b + if(x == *first2 || y == *first2){ + ++first2; + } + else{ + if ( (first1==last1) || (*first2<*first1) ) return false; + if (!(*first1<*first2)) ++first2; + ++first1; + } + } + return true; + } + + + /** + * @return true iff the simplex contains the vertex v, i.e. iff \f$ v \in (*this) \f$. + */ + bool contains(T v) const{ + return (simplex_set.find(v) != simplex_set.end()); + } + + /** + * @return \f$ (*this) \cap a = \emptyset \f$. + */ + bool disjoint(const Skeleton_blocker_simplex& a) const{ + std::vector<T> v; + v.reserve(std::min(simplex_set.size(), a.simplex_set.size())); + + set_intersection(simplex_set.cbegin(),simplex_set.cend(), + a.simplex_set.cbegin(),a.simplex_set.cend(), + std::back_inserter(v)); + + return (v.size()==0); + } + + + bool operator==(const Skeleton_blocker_simplex& other) const{ + return (this->simplex_set == other.simplex_set); + } + + bool operator!=(const Skeleton_blocker_simplex& other) const{ + return (this->simplex_set != other.simplex_set); + } + + bool operator<(const Skeleton_blocker_simplex& other) const{ + return (std::lexicographical_compare(this->simplex_set.begin(),this->simplex_set.end(), + other.begin(),other.end())); + } + + //@} + + + + + + /** + * Display a simplex + */ + friend std::ostream& operator << (std::ostream& o, const Skeleton_blocker_simplex & sigma) + { + bool first = true; + o << "{"; + for(auto i : sigma) + { + if(first) first = false ; + else o << ","; + o << i; + } + o << "}"; + return o; + } + + +}; + +} + +} // namespace GUDHI + +#endif + + diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h new file mode 100644 index 00000000..5ca8e9af --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h @@ -0,0 +1,295 @@ +#ifndef GUDHI_SKELETON_BLOCKER_SUB_COMPLEX_H +#define GUDHI_SKELETON_BLOCKER_SUB_COMPLEX_H + +#include "gudhi/Skeleton_blocker_complex.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h" +#include "gudhi/Utils.h" + +namespace Gudhi{ + +namespace skbl { + +/** + * @brief Simplicial subcomplex of a complex represented by a skeleton/blockers pair. + * + * @details Stores a subcomplex of a simplicial complex. + * To simplify explanations below, we will suppose that : + * - K is the root simplicial complex + * - L is a subcomplex of K. + * + * One vertex of K may exists in L but with a different address. + * To be able to locate the vertices in K from vertices of L, the class + * stores a map 'adresses' between vertices of K and vertices of L. + * + * Note that the type for handle of vertices of L is 'Vertex_handle' and + * the type for handle of vertices of K is 'Root_vertex_handle'. + * + * The template ComplexType is type of the root complex. It allows to know + * if the subcomplex is geometric or not. + * It has to be either 'Skeleton_blockers_complex' or 'Skeleton_blockers_geometric_complex'. + * + */ +template<typename ComplexType> +class Skeleton_blocker_sub_complex : public ComplexType +{ + +protected: + + template<class T> friend class Skeleton_blocker_link_complex; + + typedef typename ComplexType::Graph Graph; + typedef typename ComplexType::Edge_handle Edge_handle; + + typedef typename ComplexType::boost_vertex_handle boost_vertex_handle; + + +public: + using ComplexType::add_vertex; + using ComplexType::add_edge; + using ComplexType::add_blocker; + + typedef typename ComplexType::Vertex_handle Vertex_handle; + typedef typename ComplexType::Root_vertex_handle Root_vertex_handle; + typedef typename ComplexType::Simplex_handle Simplex_handle; + typedef typename ComplexType::Root_simplex_handle Root_simplex_handle; + + +protected: + + + ///** + //* @brief Returns true iff the simplex formed by all vertices contained in 'addresses_sigma_in_link' + //* but 'vertex_to_be_ignored' is in 'link' + //*/ + /* + template<typename T> friend bool + proper_face_in_union( + Skeleton_blocker_sub_complex<T> & link, + std::vector<boost::optional<typename T::Vertex_handle> > & addresses_sigma_in_link, + int vertex_to_be_ignored);*/ + + /** + * @brief Determines whether all proper faces of simplex 'sigma' belong to 'link1' \cup 'link2' + * where 'link1' and 'link2' are subcomplexes of the same complex of type ComplexType + */ + // template<typename T> friend bool + // proper_faces_in_union(Skeleton_blocker_simplex<typename T::Root_vertex_handle> & sigma, Skeleton_blocker_sub_complex<T> & link1, Skeleton_blocker_sub_complex<T> & link2){ + + //template<typename T> friend bool + //proper_faces_in_union(Skeleton_blocker_simplex<typename T::Root_vertex_handle> & sigma, Skeleton_blocker_sub_complex<T> & link1, Skeleton_blocker_sub_complex<T> & link2); + + typedef std::map<Root_vertex_handle, Vertex_handle> IdAddressMap; + typedef typename IdAddressMap::value_type AddressPair; + typedef typename IdAddressMap::iterator IdAddressMapIterator; + typedef typename IdAddressMap::const_iterator IdAddressMapConstIterator; + std::map<Root_vertex_handle, Vertex_handle> adresses; + + +public: + /** + * Add a vertex 'global' of K to L. When added to L, this vertex will receive + * another number, addresses(global), its local adress. + * return the adress where the vertex lay on L. + * The vertex corresponding to 'global' must not be already present + * in the complex. + */ + Vertex_handle add_vertex(Root_vertex_handle global){ + assert(!this->contains_vertex(global)); + Vertex_handle address(boost::add_vertex(this->skeleton)); + this->num_vertices_++; + (*this)[address].activate(); + (*this)[address].set_id(global); + adresses.insert(AddressPair(global, address)); + this->degree_.push_back(0); + return address; + } + + + /** + * Add an edge (v1_root,v2_root) to the sub-complex. + * It assumes that both vertices corresponding to v1_root and v2_root are present + * in the sub-complex. + */ + void add_edge(Root_vertex_handle v1_root, Root_vertex_handle v2_root){ + auto v1_sub(this->get_address(v1_root)); + auto v2_sub(this->get_address(v2_root)); + assert(v1_sub && v2_sub); + this->ComplexType::add_edge(*v1_sub, *v2_sub); + } + + /** + * Add a blocker to the sub-complex. + * It assumes that all vertices of blocker_root are present + * in the sub-complex. + */ + void add_blocker(const Root_simplex_handle& blocker_root){ + auto blocker_sub = this->get_address(blocker_root); + assert(blocker_sub); + this->add_blocker(new Simplex_handle(*blocker_sub)); + } + + + +public: + + /** + * Constructs the restricted complex of 'parent_complex' to + * vertices of 'simplex'. + */ + friend void make_restricted_complex(const ComplexType & parent_complex, const Simplex_handle& simplex, Skeleton_blocker_sub_complex & result){ + result.clear(); + // add vertices to the sub complex + for (auto x : simplex){ + assert(parent_complex.contains_vertex(x)); + auto x_local = result.add_vertex(parent_complex[x].get_id()); + result[x_local] = parent_complex[x]; + } + + // add edges to the sub complex + for (auto x : simplex){ + // x_neigh is the neighbor of x intersected with vertices_simplex + Simplex_handle x_neigh; + parent_complex.add_neighbours(x, x_neigh, true); + x_neigh.intersection(simplex); + for (auto y : x_neigh){ + result.add_edge(parent_complex[x].get_id(), parent_complex[y].get_id()); + } + } + + // add blockers to the sub complex + for (auto blocker : parent_complex.const_blocker_range()){ + // check if it is the first time we encounter the blocker + if (simplex.contains(*blocker)){ + Root_simplex_handle blocker_root(parent_complex.get_id(*(blocker))); + Simplex_handle blocker_restr(*result.get_simplex_address(blocker_root)); + result.add_blocker(new Simplex_handle(blocker_restr)); + } + } + } + + + + void clear(){ + adresses.clear(); + ComplexType::clear(); + } + + /** + * Compute the local vertex in L corresponding to the vertex global in K. + * runs in O(log n) if n = num_vertices() + */ + boost::optional<Vertex_handle> get_address(Root_vertex_handle global) const{ + boost::optional<Vertex_handle> res; + IdAddressMapConstIterator it = adresses.find(global); + if (it == adresses.end()) res.reset(); + else res = (*it).second; + return res; + } + + // /** + // * Allocates a simplex in L corresponding to the simplex s in K + // * with its local adresses and returns an AddressSimplex. + // */ + // boost::optional<Simplex_handle> get_address(const Root_simplex_handle & s) const; + + +//private: + /** + * same as get_address except that it will return a simplex in any case. + * The vertices that were not found are not added. + */ + // @remark should be private but problem with VS + std::vector<boost::optional<Vertex_handle> > get_addresses(const Root_simplex_handle & s) const{ + std::vector<boost::optional<Vertex_handle> > res; + for (auto i : s) + { + res.push_back(get_address(i)); + } + return res; + } + +}; + + +/** + * @remark remarque perte de temps a creer un nouveau simplexe a chaque fois + * alors qu'on pourrait utiliser a la place de 'addresses_sigma_in_link' + * un simplex avec des valeurs spéciales ComplexDS::null_vertex par exemple + * pour indiquer qu'un vertex n'appartient pas au complex + */ +template<typename ComplexType> +bool proper_face_in_union( + Skeleton_blocker_sub_complex<ComplexType> & link, + std::vector<boost::optional<typename ComplexType::Vertex_handle> > & addresses_sigma_in_link, + int vertex_to_be_ignored) +{ + // we test that all vertices of 'addresses_sigma_in_link' but 'vertex_to_be_ignored' + // are in link1 if it is the case we construct the corresponding simplex + bool vertices_sigma_are_in_link = true; + typename ComplexType::Simplex_handle sigma_in_link; + for (int i = 0; i < addresses_sigma_in_link.size(); ++i){ + if (i != vertex_to_be_ignored){ + if (!addresses_sigma_in_link[i]){ + vertices_sigma_are_in_link = false; + break; + } + else sigma_in_link.add_vertex(*addresses_sigma_in_link[i]); + } + } + // If one of vertices of the simplex is not in the complex then it returns false + // Otherwise, it tests if the simplex is in the complex + return vertices_sigma_are_in_link && link.contains(sigma_in_link); +} + +/* + template<typename ComplexType> + bool + proper_faces_in_union(Skeleton_blocker_simplex<typename ComplexType::Root_vertex_handle> & sigma, Skeleton_blocker_sub_complex<ComplexType> & link1, Skeleton_blocker_sub_complex<ComplexType> & link2) + { + typedef typename ComplexType::Vertex_handle Vertex_handle; + std::vector<boost::optional<Vertex_handle> > addresses_sigma_in_link1 = link1.get_addresses(sigma); + std::vector<boost::optional<Vertex_handle> > addresses_sigma_in_link2 = link2.get_addresses(sigma); + + for (int current_index = 0; current_index < addresses_sigma_in_link1.size(); ++current_index) + { + + if (!proper_face_in_union(link1, addresses_sigma_in_link1, current_index) + && !proper_face_in_union(link2, addresses_sigma_in_link2, current_index)){ + return false; + } + } + return true; + }*/ + + + +// Remark: this function should be friend in order to leave get_adresses private +// however doing so seemes currently not possible due to a visual studio bug c2668 +// "the compiler does not support partial ordering of template functions as specified in the C++ Standard" +// http://www.serkey.com/error-c2668-ambiguous-call-to-overloaded-function-bb45ft.html +template<typename ComplexType> +bool +proper_faces_in_union(Skeleton_blocker_simplex<typename ComplexType::Root_vertex_handle> & sigma, Skeleton_blocker_sub_complex<ComplexType> & link1, Skeleton_blocker_sub_complex<ComplexType> & link2) +{ + typedef typename ComplexType::Vertex_handle Vertex_handle; + std::vector<boost::optional<Vertex_handle> > addresses_sigma_in_link1 = link1.get_addresses(sigma); + std::vector<boost::optional<Vertex_handle> > addresses_sigma_in_link2 = link2.get_addresses(sigma); + + for (int current_index = 0; current_index < addresses_sigma_in_link1.size(); ++current_index) + { + + if (!proper_face_in_union(link1, addresses_sigma_in_link1, current_index) + && !proper_face_in_union(link2, addresses_sigma_in_link2, current_index)){ + return false; + } + } + return true; +} + +} // namespace skbl + +} // namespace GUDHI + + +#endif + diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h new file mode 100644 index 00000000..0d870f1a --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h @@ -0,0 +1,53 @@ +/* + * Top_faces.h + * + * Created on: Oct 23, 2014 + * Author: dsalinas + */ + +#ifndef TOP_FACES_H_ +#define TOP_FACES_H_ + +#include <list> +#include <vector> +#include <set> + +template<typename SimplexHandle> +std::list<SimplexHandle> subfaces(SimplexHandle top_face){ + std::list<SimplexHandle> res; + if(top_face.dimension()==-1) return res; + if(top_face.dimension()==0) { + res.push_back(top_face); + return res; + } + else{ + auto first_vertex = top_face.first_vertex(); + top_face.remove_vertex(first_vertex); + res = subfaces(top_face); + std::list<SimplexHandle> copy = res; + for(auto& simplex : copy){ + simplex.add_vertex(first_vertex); + } + res.push_back(SimplexHandle(first_vertex)); + res.splice(res.end(),copy); + return res; + } +} + +/** + * add all faces of top_face in simplices_per_dimension + */ +template<typename SimplexHandle> +void register_faces( + std::vector< std::set<SimplexHandle> >& simplices_per_dimension, + const SimplexHandle& top_face){ + std::list<SimplexHandle> subfaces_list = subfaces(top_face); + for(auto& simplex : subfaces_list ){ + simplices_per_dimension[simplex.dimension()].insert(simplex); + } +} + + + + +#endif /* TOP_FACES_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h new file mode 100644 index 00000000..21091d10 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h @@ -0,0 +1,120 @@ +/* + * Skeleton_blockers_blockers_iterators.h + * + * Created on: Aug 25, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_ +#define GUDHI_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_ + +#include "boost/iterator/iterator_facade.hpp" + +namespace Gudhi{ + +namespace skbl { + +/** + * @brief Iterator through the blockers of a vertex. + */ +// ReturnType = const Simplex_handle* or Simplex_handle* +// MapIteratorType = BlockerMapConstIterator or BlockerMapIterator +template<typename MapIteratorType, typename ReturnType> +class Blocker_iterator_internal : public boost::iterator_facade< + Blocker_iterator_internal<MapIteratorType,ReturnType>, + ReturnType, + boost::forward_traversal_tag, + ReturnType + >{ +private: + MapIteratorType current_position; + MapIteratorType end_of_map; +public: + + Blocker_iterator_internal():current_position(){} + + Blocker_iterator_internal(MapIteratorType position,MapIteratorType end_of_map_ ): + current_position(position), end_of_map(end_of_map_) + { } + + bool equal(const Blocker_iterator_internal& other) const{ + return current_position == other.current_position; + } + + void increment(){ + goto_next_blocker(); + } + + ReturnType dereference() const { + return(current_position->second); + } + +private: + /** + * Let the current pair be (v,sigma) where v is a vertex and sigma is a blocker. + * If v is not the first vertex of sigma then we already have seen sigma as a blocker + * and we look for the next one. + */ + void goto_next_blocker(){ + do { + ++current_position; + } while (!(current_position == end_of_map) && !first_time_blocker_is_seen()); + } + + bool first_time_blocker_is_seen() const{ + return current_position->first == current_position->second->first_vertex(); + } +}; + + + +/** + * @brief Iterator through the blockers of a vertex + */ +// ReturnType = const Simplex_handle* or Simplex_handle* +// MapIteratorType = BlockerMapConstIterator or BlockerMapIterator +template<typename MapIteratorType, typename ReturnType> +class Blocker_iterator_around_vertex_internal : public boost::iterator_facade< + Blocker_iterator_around_vertex_internal<MapIteratorType,ReturnType>, + ReturnType, + boost::forward_traversal_tag, + ReturnType +>{ +private: + MapIteratorType current_position_; +public: + + Blocker_iterator_around_vertex_internal():current_position_(){} + + Blocker_iterator_around_vertex_internal(MapIteratorType position): + current_position_(position) + {} + + Blocker_iterator_around_vertex_internal& operator=(Blocker_iterator_around_vertex_internal other){ + this->current_position_ = other.current_position_; + return *this; + } + + bool equal(const Blocker_iterator_around_vertex_internal& other) const{ + return current_position_ == other.current_position_; + } + + void increment(){ + current_position_++; + } + + ReturnType dereference() const{ + return(current_position_->second); + } + + + MapIteratorType current_position(){ + return this->current_position_; + } +}; + +} + +} // namespace GUDHI + +#endif /* GUDHI_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h new file mode 100644 index 00000000..82025a1c --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h @@ -0,0 +1,153 @@ +/* + * Skeleton_blockers_iterators_out.h + * + * Created on: Aug 25, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_ITERATORS_EDGES_H_ +#define GUDHI_SKELETON_BLOCKERS_ITERATORS_EDGES_H_ + +#include "boost/iterator/iterator_facade.hpp" + + +namespace Gudhi{ + +namespace skbl { + +template<typename SkeletonBlockerComplex> +class Complex_edge_around_vertex_iterator : + public boost::iterator_facade < Complex_edge_around_vertex_iterator<SkeletonBlockerComplex> + , typename SkeletonBlockerComplex::Edge_handle const + , boost::forward_traversal_tag + , typename SkeletonBlockerComplex::Edge_handle const + > +{ + friend class boost::iterator_core_access; + + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::boost_adjacency_iterator boost_adjacency_iterator; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + +private: + + const Complex* complex; + Vertex_handle v; + + boost_adjacency_iterator current_; + boost_adjacency_iterator end_; + +public: + + Complex_edge_around_vertex_iterator():complex(NULL){ + } + + Complex_edge_around_vertex_iterator(const Complex* complex_,Vertex_handle v_): + complex(complex_), + v(v_) + { + tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton); + } + + /** + * returns an iterator to the end + */ + Complex_edge_around_vertex_iterator(const Complex* complex_,Vertex_handle v_,int end): + complex(complex_), + v(v_) + { + tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton); + set_end(); + } + + bool equal(const Complex_edge_around_vertex_iterator& other) const{ + return (complex== other.complex) && (v == other.v) && (current_ == other.current_) && (end_ == other.end_); + } + + void increment(){ + if(current_ != end_) + ++current_; + } + + Edge_handle dereference() const{ + return *(*complex)[std::make_pair(v,*current_)]; + } + +private: + //remove this ugly hack + void set_end(){ + current_ = end_; + } +}; + + + +/** + *@brief Iterator on the edges of a simplicial complex. + * + */ +template<typename SkeletonBlockerComplex> +class Complex_edge_iterator : +public boost::iterator_facade < Complex_edge_iterator<SkeletonBlockerComplex> +, typename SkeletonBlockerComplex::Edge_handle const +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Edge_handle const +> + +{ + friend class boost::iterator_core_access; +public: + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::boost_edge_iterator boost_edge_iterator; + typedef typename Complex::Edge_handle Edge_handle; + + + const Complex* complex; + std::pair<boost_edge_iterator,boost_edge_iterator> edge_iterator ; + + Complex_edge_iterator():complex(NULL){ + } + + Complex_edge_iterator(const SkeletonBlockerComplex* complex_): + complex(complex_), + edge_iterator(boost::edges(complex_->skeleton)) + { + } + + /** + * return an iterator to the end + */ + Complex_edge_iterator(const SkeletonBlockerComplex* complex_,int end): + complex(complex_), + edge_iterator(boost::edges(complex_->skeleton)) + { + edge_iterator.first = edge_iterator.second; + } + + + bool equal(const Complex_edge_iterator& other) const{ + return (complex == other.complex) && (edge_iterator == other.edge_iterator); + } + + void increment(){ + if(edge_iterator.first != edge_iterator.second){ + ++(edge_iterator.first); + } + } + + Edge_handle dereference() const{ + return(*(edge_iterator.first)); + } +}; + + + +} + +} // namespace GUDHI + + +#endif /* GUDHI_SKELETON_BLOCKERS_ITERATORS_EDGES_H_ */ + + diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h new file mode 100644 index 00000000..cf827705 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h @@ -0,0 +1,21 @@ +/* + * Skeleton_blockers_iterators.h + * + * Created on: Aug 25, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_ITERATORS_H_ +#define GUDHI_SKELETON_BLOCKERS_ITERATORS_H_ + + +#include "Skeleton_blockers_vertices_iterators.h" +#include "Skeleton_blockers_edges_iterators.h" +#include "Skeleton_blockers_blockers_iterators.h" +#include "Skeleton_blockers_triangles_iterators.h" +#include "Skeleton_blockers_simplices_iterators.h" + + + + +#endif /* GUDHI_SKELETON_BLOCKERS_ITERATORS_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h new file mode 100644 index 00000000..5ac0f21e --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h @@ -0,0 +1,412 @@ +/* + * Skeleton_blockers_simplices_iterators.h + * + * Created on: Sep 26, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_KELETON_BLOCKERS_SIMPLICES_ITERATORS_H_ +#define GUDHI_SKELETON_BLOCKERS_SIMPLICES_ITERATORS_H_ + +#include <memory> +#include <list> +#include <iostream> +#include "gudhi/Utils.h" +#include "boost/iterator/iterator_facade.hpp" + + +#include "gudhi/Skeleton_blocker_link_complex.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h" + + + + +namespace Gudhi { + + +namespace skbl { + + +/** + * Link may be Skeleton_blocker_link_complex<SkeletonBlockerComplex> to iterate over all + * simplices around a vertex OR + * Skeleton_blocker_superior_link_complex<SkeletonBlockerComplex> to iterate over all + * superior vertices around a vertex. + * The iteration is done by computing a trie with the link and doing a breadth-first traversal + * of the trie. + */ +template<typename SkeletonBlockerComplex,typename Link> +class Simplex_around_vertex_iterator : + public boost::iterator_facade < Simplex_around_vertex_iterator<SkeletonBlockerComplex,Link> +, typename SkeletonBlockerComplex::Simplex_handle +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Simplex_handle +> +{ + friend class boost::iterator_core_access; + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + typedef typename Complex::Simplex_handle Simplex_handle; + + + typedef typename Link::Vertex_handle Link_vertex_handle; + // Link_vertex_handle == Complex_Vertex_handle but this renaming helps avoiding confusion + + +private: + + struct Trie{ + Vertex_handle v; + std::vector<std::shared_ptr<Trie> > childs; + private: + const Trie* parent_; + public: + + + //std::list<std::unique_ptr<Trie> > childs; -> use of deleted function + + Trie():parent_(0){} + Trie(Vertex_handle v_):v(v_),parent_(0){ + } + + Trie(Vertex_handle v_,Trie* parent):v(v_),parent_(parent){ + } + + + bool operator==(const Trie& other) const{ + return (v == other.v) ; + } + + void add_child(Trie* child){ + if(child){ + std::shared_ptr<Trie> ptr_to_add(child); + childs.push_back(ptr_to_add); + child->parent_ = this; + } + } + + + friend std::ostream& operator<<(std::ostream& stream, const Trie& trie){ + stream<< "T( "<< trie.v<< " "; + for(auto t : trie.childs) + stream << *t ; + stream<<")"; + return stream; + } + + // goes to the root in the trie to consitute simplex + void add_vertices_up_to_the_root(Simplex_handle& res) const{ + res.add_vertex(v); + if(parent_) + parent_->add_vertices_up_to_the_root(res); + } + + Simplex_handle simplex() const{ + Simplex_handle res; + add_vertices_up_to_the_root(res); + return res; + } + + bool is_leaf() const{ + return childs.empty(); + } + + bool is_root() const{ + return parent_==0; + } + + const Trie* parent() { + return parent_; + } + + void remove_leaf() { + assert(is_leaf); + if(!is_root()) + parent_->childs.erase(this); + } + + private: + + + public: + + Trie* go_bottom_left(){ + if(is_leaf()) + return this; + else + return (*childs.begin())->go_bottom_left(); + } + + }; + +private: + const Complex* complex; + Vertex_handle v; + std::shared_ptr<Link> link_v; + std::shared_ptr<Trie> trie; + std::list<Trie*> nodes_to_be_seen; // todo regrouper + +public: + Simplex_around_vertex_iterator():complex(0){ + } + + Simplex_around_vertex_iterator(const Complex* complex_,Vertex_handle v_): + complex(complex_), + v(v_), + link_v(new Link(*complex_,v_)), + trie(new Trie(v_)){ + compute_trie_and_nodes_to_be_seen(); + } + + // todo avoid useless copy + // todo currently just work if copy begin iterator + Simplex_around_vertex_iterator(const Simplex_around_vertex_iterator& other): + complex(other.complex), + v(other.v), + link_v(other.link_v), + trie(other.trie), + nodes_to_be_seen(other.nodes_to_be_seen) + { + if(!other.is_end()){ + } + } + + /** + * returns an iterator to the end + */ + Simplex_around_vertex_iterator(const Complex* complex_,Vertex_handle v_,bool end): + complex(complex_), + v(v_){ + set_end(); + } + +private: + + + void compute_trie_and_nodes_to_be_seen(){ + // once we go through every simplices passing through v0 + // we remove v0. That way, it prevents from passing a lot of times + // though edges leaving v0. + // another solution would have been to provides an adjacency iterator + // to superior vertices that avoids lower ones. + while(!link_v->empty()){ + auto v0 = *(link_v->vertex_range().begin()); + trie->add_child(build_trie(v0,trie.get())); + link_v->remove_vertex(v0); + } + nodes_to_be_seen.push_back(trie.get()); + } + + Trie* build_trie(Link_vertex_handle link_vh,Trie* parent){ + Trie* res = new Trie(parent_vertex(link_vh),parent); + for(Link_vertex_handle nv : link_v->vertex_range(link_vh)) { + if(link_vh < nv){ + Simplex_handle simplex_node_plus_nv(res->simplex()); + simplex_node_plus_nv.add_vertex(parent_vertex(nv)); + if(complex->contains(simplex_node_plus_nv)){ + res->add_child(build_trie(nv,res)); + } + } + } + return res; + } + + bool is_node_in_complex(Trie* trie){ + return true; + } + + Vertex_handle parent_vertex(Link_vertex_handle link_vh) const{ + return complex->convert_handle_from_another_complex(*link_v,link_vh); + } + + + +public: + + friend std::ostream& operator<<(std::ostream& stream, const Simplex_around_vertex_iterator& savi){ + stream<< savi.trie<< std::endl; ; + stream << "("<<savi.nodes_to_be_seen.size()<<") nodes to see\n"; + return stream; + } + + bool equal(const Simplex_around_vertex_iterator& other) const{ + bool same_complex = (complex == other.complex); + if(!same_complex) + return false; + + if(!(v == other.v)) + return false; + + bool both_empty = nodes_to_be_seen.empty() && other.nodes_to_be_seen.empty(); + if(both_empty) + return true; + + bool both_non_empty = !nodes_to_be_seen.empty() && !other.nodes_to_be_seen.empty(); + + if(!both_non_empty) return false; //one is empty the other is not + + bool same_node = (**(nodes_to_be_seen.begin()) == **(other.nodes_to_be_seen.begin())); + return same_node; + } + + void increment(){ + assert(!is_end()); + Trie* first_node = nodes_to_be_seen.front(); + + nodes_to_be_seen.pop_front(); + + for(auto childs : first_node->childs){ + nodes_to_be_seen.push_back(childs.get()); + } + + } + + Simplex_handle dereference() const{ + assert(!nodes_to_be_seen.empty()); + Trie* first_node = nodes_to_be_seen.front(); + return first_node->simplex(); + } + +private: + void set_end(){ + nodes_to_be_seen.clear(); + } + + bool is_end() const{ + return nodes_to_be_seen.empty(); + } +}; + + + +template<typename SkeletonBlockerComplex> +class Simplex_iterator : + public boost::iterator_facade < Simplex_iterator<SkeletonBlockerComplex> +, typename SkeletonBlockerComplex::Simplex_handle +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Simplex_handle +> +{ + typedef Skeleton_blocker_link_superior<SkeletonBlockerComplex> Link; + + friend class boost::iterator_core_access; + + template<class SkBlDS> friend class Skeleton_blocker_complex; + + + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + typedef typename Complex::Simplex_handle Simplex_handle; + + typedef typename Complex::CVI CVI; + + + typedef typename Link::Vertex_handle Link_vertex_handle; + +private: + + const Complex* complex_; + CVI current_vertex_; + + typedef Simplex_around_vertex_iterator<SkeletonBlockerComplex,Link> SAVI; + SAVI current_simplex_around_current_vertex_; + SAVI simplices_around_current_vertex_end_; + + +public: + Simplex_iterator():complex_(0){} + + // should not be called if the complex is empty + Simplex_iterator(const Complex* complex): + complex_(complex), + current_vertex_(complex->vertex_range().begin()), + current_simplex_around_current_vertex_(complex,*current_vertex_), + simplices_around_current_vertex_end_(complex,*current_vertex_,true) + { + assert(!complex->empty()); + } + +private: + // todo return to private + Simplex_iterator(const Complex* complex,bool end): + complex_(complex) + { + set_end(); + } + +public: + + Simplex_iterator(const Simplex_iterator& other) + : + complex_(other.complex_), + current_vertex_(other.current_vertex_), + current_simplex_around_current_vertex_(other.current_simplex_around_current_vertex_), + simplices_around_current_vertex_end_(other.simplices_around_current_vertex_end_) + { + } + + friend Simplex_iterator make_begin_iterator(const Complex* complex){ + if(complex->empty()) + return make_end_simplex_iterator(complex); + else + return Simplex_iterator(complex); + } + + friend Simplex_iterator make_end_simplex_iterator(const Complex* complex){ + return Simplex_iterator(complex,true); + } + + bool equal(const Simplex_iterator& other) const{ + if(complex_!=other.complex_) return false; + if(current_vertex_!=other.current_vertex_) return false; + if(is_end() && other.is_end()) return true; + if(current_simplex_around_current_vertex_ != other.current_simplex_around_current_vertex_) + return false; + return true; + } + + void increment(){ + if(current_simplex_around_current_vertex_!= simplices_around_current_vertex_end_){ + current_simplex_around_current_vertex_.increment(); + if( current_simplex_around_current_vertex_== simplices_around_current_vertex_end_) + goto_next_vertex(); + } + else{ + goto_next_vertex(); + } + } + + void goto_next_vertex(){ + current_vertex_.increment(); + if(!is_end()){ + current_simplex_around_current_vertex_= SAVI(complex_,*current_vertex_); + simplices_around_current_vertex_end_ = SAVI(complex_,*current_vertex_,true); + } + } + + Simplex_handle dereference() const{ + return current_simplex_around_current_vertex_.dereference(); + } + +private: + void set_end(){ + current_vertex_ = complex_->vertex_range().end(); + } + + bool is_end() const{ + return (current_vertex_ == complex_->vertex_range().end()); + } +}; + + +} + +} // namespace GUDHI + + + + + +#endif /* SKELETON_BLOCKERS_SIMPLICES_ITERATORS_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h new file mode 100644 index 00000000..4b66ced5 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h @@ -0,0 +1,224 @@ +/* + * Skeleton_blockers_triangles_iterators.h + * + * Created on: Aug 25, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_ +#define GUDHI_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_ + +#include "boost/iterator/iterator_facade.hpp" +#include <memory> + +namespace Gudhi{ + +namespace skbl { + +////////////////////////////////////////////////////////////////////// +/** + * \brief Iterator over the triangles that are + * adjacent to a vertex of the simplicial complex. + * \remark Will be removed soon -> dont look + */ +template<typename Complex,typename LinkType> +class Triangle_around_vertex_iterator : public boost::iterator_facade +< Triangle_around_vertex_iterator <Complex,LinkType> +, typename Complex::Simplex_handle const +, boost::forward_traversal_tag +, typename Complex::Simplex_handle const +> +{ + friend class boost::iterator_core_access; + template<typename T> friend class Triangle_iterator; +private: + typedef typename LinkType::Vertex_handle Vertex_handle; + typedef typename LinkType::Root_vertex_handle Root_vertex_handle; + typedef typename LinkType::Simplex_handle Simplex_handle; + typedef Complex_edge_iterator<Complex> Complex_edge_iterator_; + + const Complex* complex_; + Vertex_handle v_; + std::shared_ptr<LinkType> link_; + Complex_edge_iterator_ current_edge_; + bool is_end_; +public: + Triangle_around_vertex_iterator(const Complex* complex,Vertex_handle v): + complex_(complex),v_(v),link_(new LinkType(*complex,v_)), + current_edge_(link_->edge_range().begin()), + is_end_(current_edge_ == link_->edge_range().end()){ + } + + /** + * @brief ugly hack to get an iterator to the end + */ + Triangle_around_vertex_iterator(const Complex* complex,Vertex_handle v,bool is_end): + complex_(complex),v_(v),link_(0),is_end_(true){ + } + + /** + * @brief ugly hack to get an iterator to the end + */ + Triangle_around_vertex_iterator(): + complex_(0),v_(-1),link_(0),is_end_(true){ + } + + + Triangle_around_vertex_iterator(const Triangle_around_vertex_iterator& other){ + v_ = other.v_; + complex_ = other.complex_; + is_end_ = other.is_end_; + + if(!is_end_){ + link_ = other.link_; + current_edge_= other.current_edge_; + } + } + + bool equal(const Triangle_around_vertex_iterator& other) const{ + return (complex_==other.complex_) && ((finished() &&other.finished()) || current_edge_ == other.current_edge_); + } + + Simplex_handle dereference() const{ + Root_vertex_handle v1 = (*link_)[*current_edge_].first(); + Root_vertex_handle v2 = (*link_)[*current_edge_].second(); + return Simplex_handle(v_,*(complex_->get_address(v1)),*(complex_->get_address(v2))); + } + + void increment(){ + ++current_edge_; + } + +private: + bool finished() const{ + return is_end_ || (current_edge_ == link_->edge_range().end()); + } + +}; + + + +/** + * \brief Iterator over the triangles of the + * simplicial complex. + * \remark Will be removed soon -> dont look + * + */ +template<typename SkeletonBlockerComplex> +class Triangle_iterator : + public boost::iterator_facade< + Triangle_iterator <SkeletonBlockerComplex>, + typename SkeletonBlockerComplex::Simplex_handle const + , boost::forward_traversal_tag + , typename SkeletonBlockerComplex::Simplex_handle const + > +{ + friend class boost::iterator_core_access; +private: + typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle; + typedef typename SkeletonBlockerComplex::Root_vertex_handle Root_vertex_handle; + typedef typename SkeletonBlockerComplex::Simplex_handle Simplex_handle; + typedef typename SkeletonBlockerComplex::Superior_triangle_around_vertex_iterator STAVI; + + const SkeletonBlockerComplex* complex_; + Complex_vertex_iterator<SkeletonBlockerComplex> current_vertex_; + STAVI current_triangle_; + bool is_end_; +public: + + /* + * @remark assume that the complex is non-empty + */ + Triangle_iterator(const SkeletonBlockerComplex* complex): + complex_(complex), + current_vertex_(complex->vertex_range().begin()), + current_triangle_(complex,*current_vertex_), // xxx this line is problematic is the complex is empty + is_end_(false){ + + assert(!complex->empty()); + gotoFirstTriangle(); + } + +private: + //goto to the first triangle or to the end if none + void gotoFirstTriangle(){ + if(!is_finished() && current_triangle_.finished()){ + goto_next_vertex(); + } + } +public: + + /** + * @brief ugly hack to get an iterator to the end + * @remark assume that the complex is non-empty + */ + Triangle_iterator(const SkeletonBlockerComplex* complex,bool is_end): + complex_(complex), + current_vertex_(complex->vertex_range().end()), + current_triangle_(), // xxx this line is problematic is the complex is empty + is_end_(true){ + } + + + Triangle_iterator& operator=(const Triangle_iterator & other){ + complex_ = other.complex_; + Complex_vertex_iterator<SkeletonBlockerComplex> current_vertex_; + STAVI current_triangle_; + return *this; + } + + + bool equal(const Triangle_iterator& other) const{ + bool both_are_finished = is_finished() && other.is_finished(); + bool both_arent_finished = !is_finished() && !other.is_finished(); + // if the two iterators are not finished, they must have the same state + return (complex_==other.complex_) && + (both_are_finished || + ( (both_arent_finished) && current_vertex_ == other.current_vertex_ && current_triangle_ == other.current_triangle_)); + + } + + Simplex_handle dereference() const{ + return *current_triangle_; + } + +private: + + // goto the next vertex that has a triangle pending or the + // end vertex iterator if none exists + void goto_next_vertex(){ + assert(current_triangle_.finished()); //we mush have consume all triangles passing through the vertex + assert(!is_finished()); // we must not be done + + ++current_vertex_; + + if(!is_finished()){ + current_triangle_ = STAVI(complex_, *current_vertex_); + if(current_triangle_.finished()) + goto_next_vertex(); + } + } +public: + void increment(){ + if(!current_triangle_.finished()){ + ++current_triangle_; // problem here + if(current_triangle_.finished()) + goto_next_vertex(); + } + else{ + assert(!is_finished()); + goto_next_vertex(); + } + } + +private: + bool is_finished() const{ + return is_end_ || current_vertex_ == complex_->vertex_range().end(); + } +}; + +} + +} // namespace GUDHI + +#endif /* GUDHI_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h new file mode 100644 index 00000000..875e915a --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h @@ -0,0 +1,168 @@ +/* + * Skeleton_blockers_iterators_out.h + * + * Created on: Aug 25, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_ +#define GUDHI_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_ + +#include "boost/iterator/iterator_facade.hpp" + + +namespace Gudhi{ + +namespace skbl { + +/** + *@brief Iterator on the vertices of a simplicial complex + *@remark The increment operator go to the next active vertex. + *@remark Incrementation increases Vertex_handle. + */ +template<typename SkeletonBlockerComplex> +class Complex_vertex_iterator : public boost::iterator_facade +< Complex_vertex_iterator <SkeletonBlockerComplex> + , typename SkeletonBlockerComplex::Vertex_handle const + , boost::forward_traversal_tag + , typename SkeletonBlockerComplex::Vertex_handle const + > +{ + friend class boost::iterator_core_access; + + typedef typename SkeletonBlockerComplex::boost_vertex_iterator boost_vertex_iterator; + typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle; +private: + const SkeletonBlockerComplex* complex; + std::pair<boost_vertex_iterator,boost_vertex_iterator> vertexIterator; + + +public: + Complex_vertex_iterator():complex(NULL){ + } + + Complex_vertex_iterator(const SkeletonBlockerComplex* complex_): + complex(complex_), + vertexIterator(vertices(complex_->skeleton)){ + if(!finished() && !is_active()) { + goto_next_valid(); + } + } + + /** + * return an iterator to the end. + */ + Complex_vertex_iterator(const SkeletonBlockerComplex* complex_,int end): + complex(complex_),vertexIterator(vertices(complex_->skeleton)){ + vertexIterator.first = vertexIterator.second ; + } + +public: + void increment () {goto_next_valid();} + Vertex_handle dereference() const { + return(Vertex_handle(*(vertexIterator.first))); + } + + bool equal(const Complex_vertex_iterator& other) const{ + return vertexIterator == other.vertexIterator && complex == other.complex; + } + + bool operator<(const Complex_vertex_iterator& other) const{ + return dereference()<other.dereference(); + } + +private: + bool finished() const{ + return vertexIterator.first == vertexIterator.second; + } + + void goto_next_valid(){ + ++vertexIterator.first; + if(!finished() && !is_active()){ + goto_next_valid(); + } + } + + bool is_active() const{ + return ((*complex)[Vertex_handle(*vertexIterator.first)]).is_active(); + } + +}; + + + + +template<typename SkeletonBlockerComplex> +class Complex_neighbors_vertices_iterator +: public boost::iterator_facade < Complex_neighbors_vertices_iterator<SkeletonBlockerComplex> + , typename SkeletonBlockerComplex::Vertex_handle const + , boost::forward_traversal_tag + , typename SkeletonBlockerComplex::Vertex_handle const + > +{ + friend class boost::iterator_core_access; + + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::boost_adjacency_iterator boost_adjacency_iterator; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + +private: + + const Complex* complex; + Vertex_handle v; + + boost_adjacency_iterator current_; + boost_adjacency_iterator end_; + +public: + // boost_adjacency_iterator ai, ai_end; + // for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; ++ai){ + + Complex_neighbors_vertices_iterator():complex(NULL){ + } + + Complex_neighbors_vertices_iterator(const Complex* complex_,Vertex_handle v_): + complex(complex_), + v(v_){ + tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton); + } + + /** + * returns an iterator to the end + */ + Complex_neighbors_vertices_iterator(const Complex* complex_,Vertex_handle v_,int end): + complex(complex_), + v(v_){ + tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton); + set_end(); + } + + + void increment () { + if(current_ != end_) + ++current_; + } + + Vertex_handle dereference() const { + return(Vertex_handle(*current_)); + } + + bool equal(const Complex_neighbors_vertices_iterator& other) const{ + return (complex== other.complex) && (v == other.v) && (current_ == other.current_) && (end_ == other.end_); + } + +private: + //todo remove this ugly hack + void set_end(){ + current_ = end_; + } +}; + +} + +} // namespace GUDHI + +#endif /* GUDHI_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_ */ + + diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h new file mode 100644 index 00000000..b4c1cd9e --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h @@ -0,0 +1,1422 @@ +/* + * Skeleton_blocker_complex.h + * + * Created on: Jan, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKER_COMPLEX_H +#define GUDHI_SKELETON_BLOCKER_COMPLEX_H + + +#include <map> +#include <list> +#include <set> +#include <vector> +#include <iostream> +#include <string> +#include <fstream> +#include <sstream> +#include <memory> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <boost/range/adaptor/map.hpp> + +#include "gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h" +#include "gudhi/Skeleton_blocker_link_complex.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h" + +#include "gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h" +#include "gudhi/Skeleton_blocker/internal/Top_faces.h" +#include "gudhi/Utils.h" + + +namespace Gudhi { + +namespace skbl { + + +/** + *@class Skeleton_blocker_complex + *@brief Abstract Simplicial Complex represented with a skeleton/blockers pair. + * + * A simplicial complex is completely defined by : + * - the graph of its 1-skeleton; + * - its set of blockers. + * + * The graph is a boost graph templated with SkeletonBlockerDS::Graph_vertex and SkeletonBlockerDS::Graph_edge. + * + * One can accesses to vertices via SkeletonBlockerDS::Vertex_handle, to edges via Skeleton_blocker_complex::Edge_handle and + * simplices via Skeleton_blocker_simplex<Vertex_handle>. + * + * The SkeletonBlockerDS::Root_vertex_handle serves in the case of a subcomplex (see class Skeleton_blocker_sub_complex) + * to access to the address of one vertex in the parent complex. + * + * @todo TODO Simplex_handle are not classic handle + * + */ +template<class SkeletonBlockerDS> +class Skeleton_blocker_complex +{ + template<class ComplexType> friend class Complex_vertex_iterator; + template<class ComplexType> friend class Complex_neighbors_vertices_iterator; + template<class ComplexType> friend class Complex_edge_iterator; + template<class ComplexType> friend class Complex_edge_around_vertex_iterator; + + + template<class ComplexType> friend class Skeleton_blocker_link_complex; + template<class ComplexType> friend class Skeleton_blocker_link_superior; + template<class ComplexType> friend class Skeleton_blocker_sub_complex; + +public: + + + typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex; + typedef typename SkeletonBlockerDS::Graph_edge Graph_edge; + + typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle; + typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle; + typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle; + + + typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle; + typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle; + + + typedef Simplex_handle* Blocker_handle; + + + + typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator; + typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator; + + +protected: + + typedef typename boost::adjacency_list + < boost::setS, //edges + boost::vecS, // vertices + boost::undirectedS, + Graph_vertex, + Graph_edge + > Graph; + //todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink + // (eg Link with greater vertices) + // that burdens simplex iteration / complex initialization via list of simplices. + // to avoid that, one should modify the graph by storing two lists of adjacency for every + // vertex, the one with superior and the one with lower vertices, that way there is + // no more extra cost for computation of SuperiorLink + typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator; + typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator; + +protected: + typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator; + +public: + /** + * Handle to an edge of the complex. + */ + typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle; + + + + typedef std::multimap<Vertex_handle,Simplex_handle *> BlockerMap; + typedef typename std::multimap<Vertex_handle,Simplex_handle *>::value_type BlockerPair; + typedef typename std::multimap<Vertex_handle,Simplex_handle *>::iterator BlockerMapIterator; + typedef typename std::multimap<Vertex_handle,Simplex_handle *>::const_iterator BlockerMapConstIterator; + +protected: + int num_vertices_; + int num_blockers_; + + typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor; + // typedef Visitor* Visitor_ptr; + Visitor* visitor; + + /** + * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree. + * + * This quantity is updated when adding/removing edge. + * + * This is useful because the operation + * list.size() is done in linear time. + */ + std::vector<boost_vertex_handle> degree_; + Graph skeleton; /** 1-skeleton of the simplicial complex. */ + + + //todo remove!!! +public: + + /** Each vertex can access to the blockers passing through it. */ + BlockerMap blocker_map_; + + + + +public: + + + + + ///////////////////////////////////////////////////////////////////////////// + /** @name Constructors / Destructors / Initialization + */ + //@{ + Skeleton_blocker_complex(int num_vertices_ = 0,Visitor* visitor_=NULL):visitor(visitor_){ + clear(); + for (int i=0; i<num_vertices_; ++i){ + add_vertex(); + } + } + +private: + + /** + * this nested class is used for the next constructor that takes as an input a list of simplices. + * It stores a vector where the ith cell contains a set of i-simplices + * + */ + class Simplices_sets_from_list{ + private: + typedef std::set< Simplex_handle> Container_simplices; + typedef typename Container_simplices::iterator Simplices_iterator; + int dimension_; + std::vector<Container_simplices > simplices_; + + public: + Simplices_sets_from_list(std::list<Simplex_handle>& simplices): + dimension_(-1){ + assert(!simplices.empty()); + + for(auto simplex = simplices.begin() ; simplex != simplices.end(); ++simplex ){ + dimension_ = std::max(dimension_,(int)simplex->dimension()); + } + simplices_ = std::vector<Container_simplices >(dimension_+1); + + // compute k-simplices + for(auto simplex = simplices.begin() ; simplex != simplices.end(); ++simplex ){ + simplices_[simplex->dimension()].insert(*simplex); + } + } + + Simplices_iterator begin(int k){ + assert(0<= k && k<= dimension_); + return simplices_[k].begin(); + } + + Simplices_iterator end(int k){ + assert(0<= k && k<= dimension_); + return simplices_[k].end(); + } + + + Container_simplices& simplices(int k){ + return simplices_[k]; + } + + int dimension(){ + return dimension_; + } + + bool contains(const Simplex_handle& simplex) const{ + if(simplex.dimension()>dimension_) + return false; + else + return simplices_[simplex.dimension()].find(simplex)!= simplices_[simplex.dimension()].end(); + } + + void print(){ + for(int i = 0; i < dimension_; ++i){ + std::cout << i<<"-simplices"<<std::endl; + auto l = simplices_[i]; + for(auto s : l){ + std::cout << s<<std::endl; + } + } + } + }; + + + void compute_next_expand( + Simplices_sets_from_list& simplices, + int dim, + std::list<Simplex_handle>& next_expand) + { + next_expand.clear(); + + // high_resolution_clock::time_point tbeginfor = high_resolution_clock::now(); + // auto durationlooplink = std::chrono::duration_cast<std::chrono::microseconds>( tbeginfor - tbeginfor ).count(); + + for(auto sigma = simplices.begin(dim); sigma != simplices.end(dim); ++sigma){ + // high_resolution_clock::time_point tbeg = high_resolution_clock::now(); + Simplex_handle t(*sigma); + Skeleton_blocker_link_superior<Skeleton_blocker_complex> link(*this,t); + // xxx all time here, most likely because accessing superior edges needs passing through lower one + // currently + + // durationlooplink += std::chrono::duration_cast<std::chrono::microseconds>( high_resolution_clock::now() - tbeg ).count(); + + for(auto v : link.vertex_range()){ + Vertex_handle v_in_complex(*this->get_address( link.get_id(v)) ); + t.add_vertex(v_in_complex); + next_expand.push_back(t); + t.remove_vertex(v_in_complex); + } + } + // high_resolution_clock::time_point t2 = high_resolution_clock::now(); + // auto durationlooptotal = std::chrono::duration_cast<std::chrono::microseconds>( t2 - tbeginfor ).count(); + // DBGVALUE(durationlooptotal); + // DBGVALUE(durationlooplink); + } + + + +public: + /** + * @brief Constructor with a list of simplices + * @details The list of simplices must be the list + * of simplices of a simplicial complex. + *todo rewrite as iterators + */ + Skeleton_blocker_complex(std::list<Simplex_handle>& simplices,Visitor* visitor_=NULL): + num_vertices_(0),num_blockers_(0), + visitor(visitor_){ + Simplices_sets_from_list set_simplices(simplices); + + int dim = set_simplices.dimension(); + + + // add 1-skeleton to the complex + for(auto v_it = set_simplices.begin(0); v_it != set_simplices.end(0); ++v_it) + add_vertex(); + + for(auto e_it = set_simplices.begin(1); e_it != set_simplices.end(1); ++e_it){ + Vertex_handle a = e_it->first_vertex(); + Vertex_handle b = e_it->last_vertex(); + assert(contains_vertex(a) && contains_vertex(b)); + add_edge(a,b); + } + + // then add blockers + for(int current_dim = 1 ; current_dim <=dim ; ++current_dim){ + std::list<Simplex_handle> expansion_simplices; + compute_next_expand(set_simplices,current_dim,expansion_simplices); + + for(const auto &simplex : expansion_simplices) { + if(!set_simplices.contains(simplex)){ + add_blocker(simplex); + } + } + } + } + + + + // We cannot use the default copy constructor since we need + // to make a copy of each of the blockers + Skeleton_blocker_complex(const Skeleton_blocker_complex& copy){ + visitor = NULL; + degree_ = copy.degree_; + skeleton = Graph(copy.skeleton); + num_vertices_ = copy.num_vertices_; + + num_blockers_ = 0; + // we copy the blockers + for (auto blocker : copy.const_blocker_range()){ + add_blocker(*blocker); + } + } + + Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy){ + clear(); + visitor = NULL; + degree_ = copy.degree_; + skeleton = Graph(copy.skeleton); + num_vertices_ = copy.num_vertices_; + + num_blockers_ = 0; + // we copy the blockers + for (auto blocker : copy.const_blocker_range()) + add_blocker(*blocker); + return *this; + } + + + /** + * The destructor delete all blockers allocated. + */ + virtual ~Skeleton_blocker_complex(){ + clear(); + } + + /** + * Clears the simplicial complex. After a call to this function, + * blockers are destroyed. The 1-skeleton and the set of blockers + * are both empty. + */ + virtual void clear(){ + // xxx for now the responsabilty of freeing the visitor is for + // the user + visitor = NULL; + + degree_.clear(); + num_vertices_ =0; + + remove_blockers(); + + skeleton.clear(); + } + + void set_visitor(Visitor* other_visitor){ + visitor = other_visitor; + } + + //@} + + + + + ///////////////////////////////////////////////////////////////////////////// + /** @name Vertices operations + */ + //@{ + +public: + + /** + * @brief Return a local Vertex_handle of a vertex given a global one. + * @remark Assume that the vertex is present in the complex. + */ + Vertex_handle operator[](Root_vertex_handle global) const{ + auto local(get_address(global)); + assert(local); + return *local; + } + + Graph_vertex& operator[](Vertex_handle address){ + assert(0<=address.vertex && address.vertex< boost::num_vertices(skeleton)); + return skeleton[address.vertex]; + } + + const Graph_vertex& operator[](Vertex_handle address) const{ + assert(0<=address.vertex && address.vertex< boost::num_vertices(skeleton)); + return skeleton[address.vertex]; + } + + /** + * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle. + */ + Vertex_handle add_vertex(){ + Vertex_handle address(boost::add_vertex(skeleton)); + num_vertices_++; + (*this)[address].activate(); + // safe since we now that we are in the root complex and the field 'address' and 'id' + // are identical for every vertices + (*this)[address].set_id(Root_vertex_handle(address.vertex)); + degree_.push_back(0); + if (visitor) visitor->on_add_vertex(address); + return address; + } + + + /** + * @brief Remove a vertex from the simplicial complex + * @remark In fact, it just deactivates the vertex. + */ + void remove_vertex(Vertex_handle address){ + assert(contains_vertex(address)); + // We remove b + boost::clear_vertex(address.vertex,skeleton); + (*this)[address].deactivate(); + num_vertices_--; + degree_[address.vertex]=-1; + if (visitor) visitor->on_remove_vertex(address); + } + + /** + * @return true iff the simplicial complex contains the vertex u + */ + bool contains_vertex(Vertex_handle u) const{ + if (u.vertex<0 || u.vertex>=boost::num_vertices(skeleton)) return false; + return (*this)[u].is_active(); + } + + /** + * @return true iff the simplicial complex contains the vertex u + */ + bool contains_vertex(Root_vertex_handle u) const{ + boost::optional<Vertex_handle> address = get_address(u); + return address && (*this)[*address].is_active(); + } + + + /** + * @return true iff the simplicial complex contains all vertices + * of simplex sigma + */ + bool contains_vertices(const Simplex_handle & sigma) const{ + for (auto vertex : sigma) + if(!contains_vertex(vertex)) return false; + return true; + } + + /** + * Given an Id return the address of the vertex having this Id in the complex. + * For a simplicial complex, the address is the id but it may not be the case for a SubComplex. + * + */ + virtual boost::optional<Vertex_handle> get_address(Root_vertex_handle id) const{ + boost::optional<Vertex_handle> res; + if ( id.vertex< boost::num_vertices(skeleton) ) res = Vertex_handle(id.vertex);//xxx + return res; + } + + + + /** + * return the id of a vertex of adress local present in the graph + */ + Root_vertex_handle get_id(Vertex_handle local) const{ + assert(0<=local.vertex && local.vertex< boost::num_vertices(skeleton)); + return (*this)[local].get_id(); + } + + + /** + * If the current complex is a sub (or sup) complex of 'other', it converts + * the address of a vertex v expressed in 'other' to the address of the vertex + * v in the current one. + * @remark this methods uses Root_vertex_handle to identify the vertex + */ + Vertex_handle convert_handle_from_another_complex( + const Skeleton_blocker_complex& other,Vertex_handle vh_in_other) const{ + auto vh_in_current_complex = get_address(other.get_id(vh_in_other)); + assert(vh_in_current_complex); + return *vh_in_current_complex; + } + + + int degree(Vertex_handle local) const{ + assert(0<=local.vertex && local.vertex< boost::num_vertices(skeleton)); + return degree_[local.vertex]; + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Edges operations + */ + //@{ + +public: + + /** + * @brief return an edge handle if the two vertices forms + * an edge in the complex + */ + boost::optional<Edge_handle> operator[](const std::pair<Vertex_handle,Vertex_handle>& ab) const{ + boost::optional<Edge_handle> res; + std::pair<Edge_handle,bool> edge_pair(boost::edge(ab.first.vertex,ab.second.vertex,skeleton)); + if (edge_pair.second) + res = edge_pair.first; + return res; + } + + Graph_edge& operator[](Edge_handle edge_handle){ + return skeleton[edge_handle]; + } + + const Graph_edge& operator[](Edge_handle edge_handle) const{ + return skeleton[edge_handle]; + } + + Vertex_handle first_vertex(Edge_handle edge_handle) const{ + return source(edge_handle,skeleton); + } + + Vertex_handle second_vertex(Edge_handle edge_handle) const{ + return target(edge_handle,skeleton); + } + + /** + * @brief returns the simplex made with the two vertices of the edge + */ + Simplex_handle get_vertices(Edge_handle edge_handle) const{ + auto edge((*this)[edge_handle]); + return Simplex_handle((*this)[edge.first()],(*this)[edge.second()]); + } + + /** + * @brief Adds an edge between vertices a and b + */ + Edge_handle add_edge(Vertex_handle a, Vertex_handle b){ + assert(contains_vertex(a) && contains_vertex(b)); + assert(a!=b); + + auto edge_handle((*this)[std::make_pair(a,b)]); + // std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)]; + // Edge_handle edge_descr; + // bool edge_present = pair_descr_bool.second; + if (!edge_handle) + { + edge_handle = boost::add_edge(a.vertex,b.vertex,skeleton).first; + (*this)[*edge_handle].setId(get_id(a),get_id(b)); + degree_[a.vertex]++; + degree_[b.vertex]++; + if (visitor) visitor->on_add_edge(a,b); + } + return *edge_handle; + } + + /** + * @brief Adds all edges of simplex sigma to the simplicial complex. + */ + void add_edges(const Simplex_handle & sigma){ + Simplex_handle_iterator i, j; + for (i = sigma.begin() ; i != sigma.end() ; ++i) + for (j = i, j++ ; j != sigma.end() ; ++j) + add_edge(*i,*j); + } + + /** + * @brief Removes edge ab from the simplicial complex. + */ + virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b){ + bool found; + Edge_handle edge; + tie(edge,found) = boost::edge(a.vertex,b.vertex,skeleton); + if (found) + { + if (visitor) visitor->on_remove_edge(a,b); + // if (heapCollapse.Contains(edge)) heapCollapse.Delete(edge); + boost::remove_edge(a.vertex,b.vertex,skeleton); + degree_[a.vertex]--; + degree_[b.vertex]--; + } + return edge; + } + + + /** + * @brief Removes edge from the simplicial complex. + */ + void remove_edge(Edge_handle edge){ + assert(contains_vertex(first_vertex(edge))); + assert(contains_vertex(second_vertex(edge))); + remove_edge(first_vertex(edge),second_vertex(edge)); + } + + + /** + * @brief The complex is reduced to its set of vertices. + * All the edges and blockers are removed. + */ + void keep_only_vertices(){ + remove_blockers(); + + for(auto u : vertex_range()){ + while (this->degree(u)> 0) + { + Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first)); + this->remove_edge(u,v); + } + } + } + + + /** + * @return true iff the simplicial complex contains an edge between + * vertices a and b + */ + bool contains_edge(Vertex_handle a, Vertex_handle b) const{ + //if (a.vertex<0 || b.vertex <0) return false; + return boost::edge(a.vertex,b.vertex,skeleton).second; + } + + + /** + * @return true iff the simplicial complex contains all vertices + * and all edges of simplex sigma + */ + bool contains_edges(const Simplex_handle & sigma) const{ + for (auto i = sigma.begin() ; i != sigma.end() ; ++i){ + if(!contains_vertex(*i)) return false; + for (auto j=i; ++j != sigma.end() ; ){ + if (!contains_edge(*i,*j)) + return false; + } + } + return true; + } + //@} + + + + ///////////////////////////////////////////////////////////////////////////// + /** @name Blockers operations + */ + //@{ + /** + * Adds the 2-blocker abc + */ + void add_blocker(Vertex_handle a, Vertex_handle b, Vertex_handle c){ + add_blocker(Simplex_handle(a,b,c)); + } + + /** + * Adds the 3-blocker abcd + */ + void add_blocker(Vertex_handle a, Vertex_handle b, Vertex_handle c, Vertex_handle d){ + add_blocker(Simplex_handle(a,b,c,d)); + } + + /** + * Adds the simplex blocker_pt to the set of blockers and + * returns a Blocker_handle toward it if was not present before. + */ + Blocker_handle add_blocker(const Simplex_handle& blocker){ + if (contains_blocker(blocker)) + { + //std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; + return 0; + } + else{ + if (visitor) visitor->on_add_blocker(blocker); + Blocker_handle blocker_pt = new Simplex_handle(blocker); + num_blockers_++; + auto vertex = blocker_pt->begin(); + while(vertex != blocker_pt->end()) + { + blocker_map_.insert(BlockerPair(*vertex,blocker_pt)); + ++vertex; + } + return blocker_pt; + } + } + +protected: + /** + * Adds the simplex s to the set of blockers + */ + void add_blocker(Blocker_handle blocker){ + if (contains_blocker(*blocker)) + { + //std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; + return; + } + else{ + if (visitor) visitor->on_add_blocker(*blocker); + num_blockers_++; + auto vertex = blocker->begin(); + while(vertex != blocker->end()) + { + blocker_map_.insert(BlockerPair(*vertex,blocker)); + ++vertex; + } + } + } + + +protected: + /** + * Removes sigma from the blocker map of vertex v + */ + void remove_blocker(const Blocker_handle sigma, Vertex_handle v){ + Complex_blocker_around_vertex_iterator blocker; + for (blocker = blocker_range(v).begin(); + blocker != blocker_range(v).end(); + ++blocker + ){ + if (*blocker == sigma) break; + } + if (*blocker != sigma){ + std::cerr << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n"; + assert(false); + } + else{ + blocker_map_.erase(blocker.current_position()); + } + } + +public: + /** + * Removes the simplex sigma from the set of blockers. + * sigma has to belongs to the set of blockers + */ + void remove_blocker(const Blocker_handle sigma){ + for (auto vertex : *sigma){ + remove_blocker(sigma,vertex); + } + num_blockers_--; + } + + + + + /** + * Remove all blockers, in other words, it expand the simplicial + * complex to the smallest flag complex that contains it. + */ + void remove_blockers(){ + // Desallocate the blockers + while (!blocker_map_.empty()){ + delete_blocker(blocker_map_.begin()->second); + } + num_blockers_ = 0; + blocker_map_.clear(); + } + +protected: + /** + * Removes the simplex sigma from the set of blockers. + * sigma has to belongs to the set of blockers + * + * @remark contrarily to delete_blockers does not call the destructor + */ + void remove_blocker(const Simplex_handle& sigma){ + assert(contains_blocker(sigma)); + for (auto vertex : sigma) + remove_blocker(sigma,vertex); + num_blockers_--; + } + + + + +public: + /** + * Removes the simplex s from the set of blockers + * and desallocate s. + */ + void delete_blocker(Blocker_handle sigma){ + if (visitor) visitor->on_delete_blocker(sigma); + remove_blocker(sigma); + delete sigma; + } + + + /** + * @return true iff s is a blocker of the simplicial complex + */ + bool contains_blocker(const Blocker_handle s) const{ + if (s->dimension()<2) + return false; + + Vertex_handle a = s->first_vertex(); + + for (auto blocker : const_blocker_range(a)){ + if ( s == *blocker ) + return true; + } + return false; + } + + /** + * @return true iff s is a blocker of the simplicial complex + */ + bool contains_blocker(const Simplex_handle & s) const{ + if (s.dimension()<2) + return false; + + Vertex_handle a = s.first_vertex(); + + for (auto blocker : const_blocker_range(a)){ + if ( s == *blocker ) + return true; + } + return false; + } + + +private: + /** + * @return true iff a blocker of the simplicial complex + * is a face of sigma. + */ + bool blocks(const Simplex_handle & sigma) const{ + + for(auto blocker : const_blocker_range()){ + if ( sigma.contains(*blocker) ) + return true; + } + return false; + } + + //@} + + + + + ///////////////////////////////////////////////////////////////////////////// + /** @name Neighbourhood access + */ + //@{ + +public: + /** + * @brief Adds to simplex n the neighbours of v: + * \f$ n \leftarrow n \cup N(v) \f$. + * + * If keep_only_superior is true then only vertices that are greater than v are added. + */ + virtual void add_neighbours(Vertex_handle v, Simplex_handle & n,bool keep_only_superior=false) const{ + boost_adjacency_iterator ai, ai_end; + for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; ++ai){ + if (keep_only_superior){ + if (*ai>v.vertex) + n.add_vertex(Vertex_handle(*ai)); + } + else + n.add_vertex(Vertex_handle(*ai)); + } + } + + /** + * @brief Add to simplex res all vertices which are + * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$. + * + * If 'keep_only_superior' is true then only vertices that are greater than alpha are added. + * + * todo revoir + * + */ + virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res,bool keep_only_superior=false) const{ + res.clear(); + // ---------------------------- + // Compute vertices in the link + // we compute the intersection of N(alpha_i) and store it in n + // ---------------------------- + auto alpha_vertex = alpha.begin(); + add_neighbours(*alpha_vertex,res,keep_only_superior); + for (alpha_vertex = (alpha.begin())++ ; alpha_vertex != alpha.end() ; ++alpha_vertex) + { + keep_neighbours(*alpha_vertex,res,keep_only_superior); + } + } + + /** + * @brief Eliminates from simplex n all vertices which are + * not neighbours of v: \f$ res \leftarrow res \cap N(v) \f$. + * + * If 'keep_only_superior' is true then only vertices that are greater than v are keeped. + * + */ + virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res,bool keep_only_superior=false) const{ + Simplex_handle nv; + add_neighbours(v,nv,keep_only_superior); + res.intersection(nv); + } + + /** + * @brief Eliminates from simplex n all vertices which are + * neighbours of v: \f$ res \leftarrow res \setminus N(v) \f$. + * + * If 'keep_only_superior' is true then only vertices that are greater than v are added. + * + */ + virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res,bool keep_only_superior=false) const{ + Simplex_handle nv; + add_neighbours(v,nv,keep_only_superior); + res.difference(nv); + } + //@} + + + ///////////////////////////////////////////////////////////////////////////// + /** @name Operations on the simplicial complex + */ + //@{ +public: + + /** + * @brief Compute the local vertices of 's' in the current complex + * If one of them is not present in the complex then the return value is uninitialized. + * + * xxx rename get_address et place un using dans sub_complex + */ + boost::optional<Simplex_handle> get_simplex_address(const Root_simplex_handle& s) const + { + boost::optional<Simplex_handle> res; + + Simplex_handle s_address; + //Root_simplex_const_iterator i; + for (auto i = s.begin() ; i != s.end() ; ++i) + { + boost::optional<Vertex_handle> address = get_address(*i); + if (!address) + return res; + else + s_address.add_vertex(*address); + } + res = s_address; + return res; + } + + /** + * @brief returns a simplex with vertices which are the id of vertices of the + * argument. + */ + Root_simplex_handle get_id(const Simplex_handle& local_simplex) const{ + Root_simplex_handle global_simplex; + for (auto x = local_simplex.begin(); x!= local_simplex.end();++x){ + global_simplex.add_vertex(get_id(*x)); + + } + return global_simplex; + + } + + + /** + * @brief returns true iff the simplex s belongs to the simplicial + * complex. + */ + virtual bool contains(const Simplex_handle & s) const{ + if (s.dimension() == -1 ) return false; + else + if (s.dimension() ==0 ){ + return contains_vertex(s.first_vertex()); + } + else + return ( contains_edges(s) && !blocks(s) ); + } + + /* + * @brief returnrs true iff the complex is empty. + */ + bool empty() const{ + return num_vertices()==0; + } + + /* + * @brief returns the number of vertices in the complex. + */ + int num_vertices() const{ + //remark boost::num_vertices(skeleton) counts deactivated vertices + return num_vertices_; + } + + /* + * @brief returns the number of edges in the complex. + * todo in O(n), cache the value + */ + int num_edges() const{ + return boost::num_edges(skeleton); + } + + /* + * @brief returns the number of blockers in the complex. + */ + int num_blockers() const{ + return num_blockers_; + } + + /* + * @brief returns true iff the graph of the 1-skeleton of the complex is complete. + */ + bool complete() const{ + return (num_vertices()*(num_vertices()-1))/2 == num_edges(); + } + + /** + * @brief returns the number of connected components in the graph of the 1-skeleton. + */ + int num_connected_components() const{ + int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices(); + std::vector<int> component(skeleton.vertex_set().size()); + return boost::connected_components(this->skeleton,&component[0]) - num_vert_collapsed; + } + + + //todo remove + // do + void keep_only_largest_cc(){ + std::vector<unsigned> component(skeleton.vertex_set().size()); + boost::connected_components(this->skeleton,&component[0]); + auto maxCC = min_element(component.begin(),component.end()); + for (unsigned i = 0; i != component.size(); ++i){ + if(component[i]!=*maxCC){ + if(this->contains_vertex(Vertex_handle(i))) + this->remove_vertex(Vertex_handle(i)); + } + } + } + + + /** + * @brief %Test if the complex is a cone. + * @details Runs in O(n) where n is the number of vertices. + */ + bool is_cone() const{ + if (num_vertices()==0) return false; + if (num_vertices()==1) return true; + for(auto vi : vertex_range()){ + //xxx todo faire une methode bool is_in_blocker(Vertex_handle) + if (blocker_map_.find(vi)==blocker_map_.end()){ + // no blocker passes through the vertex, we just need to + // check if the current vertex is linked to all others vertices of the complex + if (degree_[vi.vertex] == num_vertices()-1) + return true; + } + } + return false; + } + + //@} + + + ///////////////////////////////////////////////////////////////////////////// + /** @name Vertex iterators + */ + //@{ + + typedef Complex_vertex_iterator<Skeleton_blocker_complex> CVI; //todo rename + + // /** + // * @brief Range over the vertices of the simplicial complex. + // * Methods .begin() and .end() return a Complex_vertex_iterator. + // */ + typedef boost::iterator_range < Complex_vertex_iterator<Skeleton_blocker_complex> > Complex_vertex_range; + + /** + * @brief Returns a Complex_vertex_range over all vertices of the complex + */ + Complex_vertex_range vertex_range() const + { + auto begin = Complex_vertex_iterator<Skeleton_blocker_complex>(this); + auto end = Complex_vertex_iterator<Skeleton_blocker_complex>(this,0); + return Complex_vertex_range(begin,end); + } + + typedef boost::iterator_range < Complex_neighbors_vertices_iterator<Skeleton_blocker_complex> > Complex_neighbors_vertices_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v + */ + Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const + { + auto begin = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(this,v); + auto end = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(this,v,0); + return Complex_neighbors_vertices_range(begin,end); + } + + //@} + + + /** @name Edge iterators + */ + //@{ + + + typedef boost::iterator_range <Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> + Complex_edge_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex + */ + Complex_edge_range edge_range() const + { + auto begin = Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this); + auto end = Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,0); + return Complex_edge_range(begin,end); + } + + + + typedef boost::iterator_range <Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> + Complex_edge_around_vertex_range; + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes + * through 'v' + */ + Complex_edge_around_vertex_range edge_range(Vertex_handle v) const + { + auto begin = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,v); + auto end = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,v,0); + return Complex_edge_around_vertex_range(begin,end); + } + + + + //@} + + + /** @name Triangles iterators + */ + //@{ +private: + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link; + typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link; +public: + typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex,Superior_link> Superior_triangle_around_vertex_iterator; + + typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link> > Complex_triangle_around_vertex_range; + + /** + * @brief Range over triangles around a vertex of the simplicial complex. + * Methods .begin() and .end() return a Triangle_around_vertex_iterator. + * + */ + Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const + { + auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link>(this,v); + auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link>(this,v,0); + return Complex_triangle_around_vertex_range(begin,end); + } + + + typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range; + + /** + * @brief Range over triangles of the simplicial complex. + * Methods .begin() and .end() return a Triangle_around_vertex_iterator. + * + */ + Complex_triangle_range triangle_range() const + { + auto end = Triangle_iterator<Skeleton_blocker_complex>(this,0); + if(empty()){ + return Complex_triangle_range(end,end); + } + else{ + auto begin = Triangle_iterator<Skeleton_blocker_complex>(this); + return Complex_triangle_range(begin,end); + } + + } + + + //@} + + + + /** @name Simplices iterators + */ + //@{ + typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex,Link> Complex_simplex_around_vertex_iterator; + + /** + * @brief Range over the simplices of the simplicial complex around a vertex. + * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator. + */ + typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range; + + /** + * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex + */ + Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const + { + assert(contains_vertex(v)); + return Complex_simplex_around_vertex_range( + Complex_simplex_around_vertex_iterator(this,v), + Complex_simplex_around_vertex_iterator(this,v,true) + ); + } + + // typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator; + typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator; + + typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range; + + Complex_simplex_range simplex_range() const + { + Complex_simplex_iterator end(this,true); + if(empty()){ + return Complex_simplex_range(end,end); + } + else{ + Complex_simplex_iterator begin(this); + return Complex_simplex_range(begin,end); + } + } + + + //@} + + + /** @name Blockers iterators + */ + //@{ +private: + /** + * @brief Iterator over the blockers adjacent to a vertex + */ + typedef Blocker_iterator_around_vertex_internal< + typename std::multimap<Vertex_handle,Simplex_handle *>::iterator, + Blocker_handle> + Complex_blocker_around_vertex_iterator; + + /** + * @brief Iterator over (constant) blockers adjacent to a vertex + */ + typedef Blocker_iterator_around_vertex_internal< + typename std::multimap<Vertex_handle,Simplex_handle *>::const_iterator, + const Blocker_handle> + Const_complex_blocker_around_vertex_iterator; + + typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range; + typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range; + +public: + + + + Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) + { + auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); + auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); + return Complex_blocker_around_vertex_range(begin,end); + } + + /** + * @brief Returns a Complex_blocker_around_vertex_range over all blockers of the complex adjacent to the vertex v. + */ + Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const + { + auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); + auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); + return Const_complex_blocker_around_vertex_range(begin,end); + } + + + + +private: + + /** + * @brief Iterator over the blockers. + */ + typedef Blocker_iterator_internal< + typename std::multimap<Vertex_handle,Simplex_handle *>::iterator, + Blocker_handle> + Complex_blocker_iterator; + + /** + * @brief Iterator over the (constant) blockers. + */ + typedef Blocker_iterator_internal< + typename std::multimap<Vertex_handle,Simplex_handle *>::const_iterator, + const Blocker_handle> + Const_complex_blocker_iterator; + + typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range; + typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range; + + +public: + + + Complex_blocker_range blocker_range() + { + auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end() ); + auto end = Complex_blocker_iterator(blocker_map_.end() , blocker_map_.end() ); + return Complex_blocker_range(begin,end); + } + + + Const_complex_blocker_range const_blocker_range() const + { + auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end() ); + auto end = Const_complex_blocker_iterator(blocker_map_.end() , blocker_map_.end() ); + return Const_complex_blocker_range(begin,end); + } + + + + + //@} + + + + + + ///////////////////////////////////////////////////////////////////////////// + /** @name Print and IO methods + */ + //@{ +public: + std::string to_string() const{ + std::ostringstream stream; + stream<<num_vertices()<<" vertices:\n"<<vertices_to_string()<<std::endl; + stream<<num_edges()<<" edges:\n"<<edges_to_string()<<std::endl; + stream<<num_blockers()<<" blockers:\n"<<blockers_to_string()<<std::endl; + return stream.str(); + } + + std::string vertices_to_string() const{ + std::ostringstream stream; + for(auto vertex : vertex_range()) + stream << "("<<(*this)[vertex].get_id()<<"),"; + stream<< std::endl; + return stream.str(); + } + + std::string edges_to_string() const{ + std::ostringstream stream; + for(auto edge : edge_range()) + stream << "("<< (*this)[edge].first()<<","<< (*this)[edge].second() << ")"<<" id = "<< (*this)[edge].index()<< std::endl; + stream<< std::endl; + return stream.str(); + } + + + std::string blockers_to_string() const{ + std::ostringstream stream; + for (auto bl:blocker_map_){ + stream << bl.first << " => " << bl.second << ":"<<*bl.second <<"\n"; + } + return stream.str(); + } + + //@} + +}; + + + +/** + * build a simplicial complex from a collection + * of top faces. + */ +template<typename Complex,typename SimplexHandleIterator> +unsigned make_complex_from_top_faces(Complex& complex,SimplexHandleIterator begin,SimplexHandleIterator end){ + typedef typename Complex::Simplex_handle Simplex_handle; + + int dimension = 0; + for(auto top_face = begin; top_face != end; ++top_face) + dimension = std::max(dimension,top_face->dimension()); + + std::vector< std::set<Simplex_handle> > simplices_per_dimension(dimension+1); + + + for(auto top_face = begin; top_face != end; ++top_face){ + register_faces(simplices_per_dimension,*top_face); + } + + // compute list of simplices + std::list<Simplex_handle> simplices; + for(int dim = 0 ; dim <= dimension ; ++dim){ + std::copy( + simplices_per_dimension[dim].begin(), + simplices_per_dimension[dim].end(), + std::back_inserter(simplices) + ); + } + complex = Complex(simplices); + return simplices.size(); +} + + + +} + +} // namespace GUDHI + + + +#endif + diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h new file mode 100644 index 00000000..924b653c --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h @@ -0,0 +1,118 @@ +/* + * Skeleton_blocker_geometric_complex.h + * + * Created on: Feb 11, 2014 + * Author: dsalinas + */ + +#ifndef SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_ +#define SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_ + + +#include "gudhi/Utils.h" +#include "gudhi/Skeleton_blocker_simplifiable_complex.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h" + + + +namespace Gudhi{ + +namespace skbl { + +/** + * @brief Class that represents a geometric complex that can be simplified. + * The class allows access to points of vertices. + * + */ +template<typename SkeletonBlockerGeometricDS> +class Skeleton_blocker_geometric_complex : public Skeleton_blocker_simplifiable_complex<SkeletonBlockerGeometricDS> +{ +public: + typedef typename SkeletonBlockerGeometricDS::GT GT; + + typedef Skeleton_blocker_simplifiable_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker; + + typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle; + typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle; + typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle; + typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle; + + + typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex; + + typedef typename SkeletonBlockerGeometricDS::Point Point; + + + /** + * @brief Add a vertex to the complex with its associated point. + */ + void add_vertex(const Point& point){ + Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex(); + (*this)[ad].point() = point; + } + + + /** + * @brief Returns the Point associated to the vertex v. + */ + const Point& point(Vertex_handle v) const{ + assert(this->contains_vertex(v)); + return (*this)[v].point(); + } + + /** + * @brief Returns the Point associated to the vertex v. + */ + Point& point(Vertex_handle v) { + assert(this->contains_vertex(v)); + return (*this)[v].point(); + } + + const Point& point(Root_vertex_handle global_v) const{ + Vertex_handle local_v ( (*this)[global_v]) ; + assert(this->contains_vertex(local_v)); + return (*this)[local_v].point(); + } + + Point& point(Root_vertex_handle global_v) { + Vertex_handle local_v ( (*this)[global_v]) ; + assert(this->contains_vertex(local_v)); + return (*this)[local_v].point(); + } + + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link; + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(const Simplex_handle& simplex) const{ + Geometric_link link(*this,simplex); + //we now add the point info + add_points_to_link(link); + return link; + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(Edge_handle edge) const{ + Geometric_link link(*this,edge); + //we now add the point info + add_points_to_link(link); + return link; + } +private: + void add_points_to_link(Geometric_link& link) const{ + for(Vertex_handle v : link.vertex_range()){ + Root_vertex_handle v_root(link.get_id(v)); + link.point(v) = (*this).point(v_root); + } + } +}; + +} + +} // namespace GUDHI + +#endif /* SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h new file mode 100644 index 00000000..29c0691b --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h @@ -0,0 +1,275 @@ +/* + * Skeleton_blockers_link_complex.h + * + * Created on: Feb 7, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_LINK_COMPLEX_H_ +#define GUDHI_SKELETON_BLOCKERS_LINK_COMPLEX_H_ + +#include "gudhi/Utils.h" +#include "gudhi/Skeleton_blocker_complex.h" + +namespace Gudhi{ + +namespace skbl { + +template<class ComplexType> class Skeleton_blocker_sub_complex; + + +/** + * \brief Class representing the link of a simplicial complex encoded by a skeleton/blockers pair. + * It inherits from Skeleton_blocker_sub_complex because such complex is a sub complex of a + * root complex. + */ +template<typename ComplexType> +class Skeleton_blocker_link_complex : public Skeleton_blocker_sub_complex<ComplexType> +{ + template <typename T> friend class Skeleton_blocker_link_superior; + typedef typename ComplexType::Edge_handle Edge_handle; + + typedef typename ComplexType::boost_vertex_handle boost_vertex_handle; + + +private: + + bool only_superior_vertices_; + +public: + typedef typename ComplexType::Vertex_handle Vertex_handle; + typedef typename ComplexType::Root_vertex_handle Root_vertex_handle; + + typedef typename ComplexType::Simplex_handle Simplex_handle; + typedef typename ComplexType::Root_simplex_handle Root_simplex_handle; + + typedef typename ComplexType::Blocker_handle Blocker_handle; + + + typedef typename ComplexType::Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator; + typedef typename ComplexType::Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_handle_iterator; + + + Skeleton_blocker_link_complex(bool only_superior_vertices=false):only_superior_vertices_(only_superior_vertices){ + } + + /** + * If the parameter only_superior_vertices is true, + * only vertices greater than the one of alpha are added. + */ + Skeleton_blocker_link_complex(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool only_superior_vertices = false) + :only_superior_vertices_(only_superior_vertices) { + if(!alpha_parent_adress.empty()) + build_link(parent_complex,alpha_parent_adress); + + } + + /** + * If the parameter only_superior_vertices is true, + * only vertices greater than the one of the vertex are added. + */ + Skeleton_blocker_link_complex(const ComplexType & parent_complex, Vertex_handle a_parent_adress, bool only_superior_vertices = false) + :only_superior_vertices_(only_superior_vertices){ + Simplex_handle alpha_simplex(a_parent_adress); + build_link(parent_complex,alpha_simplex); + } + + + /** + * If the parameter only_superior_vertices is true, + * only vertices greater than the one of the edge are added. + */ + Skeleton_blocker_link_complex(const ComplexType & parent_complex, Edge_handle edge, bool only_superior_vertices = false) + :only_superior_vertices_(only_superior_vertices){ + Simplex_handle alpha_simplex(parent_complex.first_vertex(edge),parent_complex.second_vertex(edge)); + build_link(parent_complex,alpha_simplex); + } + +protected: + + + /** + * @brief compute vertices of the link. + * If the boolean only_superior_vertices is true, then only the vertices + * are greater than vertices of alpha_parent_adress are added. + */ + void compute_link_vertices(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool only_superior_vertices,bool is_alpha_blocker = false){ + if(alpha_parent_adress.dimension()==0) + // for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector) + // thus we call a specific function that will reserve a vector with appropriate size + this->compute_link_vertices(parent_complex,alpha_parent_adress.first_vertex(),only_superior_vertices_); + else{ + // we compute the intersection of neighbors of alpha and store it in link_vertices + Simplex_handle link_vertices_parent; + parent_complex.add_neighbours(alpha_parent_adress,link_vertices_parent,only_superior_vertices); + // For all vertex 'v' in this intersection, we go through all its adjacent blockers. + // If one blocker minus 'v' is included in alpha then the vertex is not in the link complex. + for (auto v_parent : link_vertices_parent){ + bool new_vertex = true; + for (auto beta : parent_complex.const_blocker_range(v_parent)) + { + if(!is_alpha_blocker || *beta!=alpha_parent_adress){ + new_vertex = !(alpha_parent_adress.contains_difference(*beta,v_parent)); + if(!new_vertex) break; + } + } + if (new_vertex) + this->add_vertex(parent_complex.get_id(v_parent)); + } + } + } + + + /** + * @brief compute vertices of the link. + * If the boolean only_superior_vertices is true, then only the vertices + * are greater than vertices of alpha_parent_adress are added. + */ + void compute_link_vertices(const ComplexType & parent_complex, Vertex_handle alpha_parent_adress,bool only_superior_vertices){ + // for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector + this->skeleton.m_vertices.reserve(parent_complex.degree(alpha_parent_adress)); + + // For all vertex 'v' in this intersection, we go through all its adjacent blockers. + // If one blocker minus 'v' is included in alpha then the vertex is not in the link complex. + for (auto v_parent : parent_complex.vertex_range(alpha_parent_adress)){ + if(!only_superior_vertices || v_parent.vertex> alpha_parent_adress.vertex) + this->add_vertex(parent_complex.get_id(v_parent)); + } + + } + + + void compute_link_edges(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool is_alpha_blocker = false){ + Simplex_handle_iterator y_link, x_parent , y_parent; + // ---------------------------- + // Compute edges in the link + // ------------------------- + + if(this->num_vertices()<=1) return; + + for (auto x_link = this->vertex_range().begin() ; + x_link!= this->vertex_range().end(); + ++x_link + ) + { + for (auto y_link = x_link; ++y_link != this->vertex_range().end(); ){ + Vertex_handle x_parent = *parent_complex.get_address(this->get_id(*x_link)); + Vertex_handle y_parent = *parent_complex.get_address(this->get_id(*y_link)); + if (parent_complex.contains_edge(x_parent,y_parent) ){ + // we check that there is no blocker subset of alpha passing trough x and y + bool new_edge=true; + for (auto blocker_parent : parent_complex.const_blocker_range(x_parent)) + { + if(!is_alpha_blocker || *blocker_parent!=alpha_parent_adress){ + if (blocker_parent->contains(y_parent)) + { + new_edge = !(alpha_parent_adress.contains_difference(*blocker_parent,x_parent,y_parent)); + if (!new_edge) break; + } + } + } + if (new_edge) + this->add_edge(*x_link,*y_link); + } + } + } + } + + + /** + * @brief : Given an address in the current complex, it returns the + * corresponding address in 'other_complex'. + * It assumes that other_complex have a vertex 'this.get_id(address)' + */ + boost::optional<Vertex_handle> give_equivalent_vertex(const ComplexType & other_complex, + Vertex_handle address) const{ + Root_vertex_handle id((*this)[address].get_id()); + return other_complex.get_address(id); + } + + /* + * compute the blockers of the link if is_alpha_blocker is false. + * Otherwise, alpha is a blocker, and the link is computed in the complex where + * the blocker is anticollapsed. + */ + void compute_link_blockers(const ComplexType & parent_complex,const Simplex_handle& alpha_parent,bool is_alpha_blocker = false){ + + for (auto x_link : this->vertex_range()){ + + Vertex_handle x_parent = * this->give_equivalent_vertex(parent_complex,x_link); + + for (auto blocker_parent : parent_complex.const_blocker_range(x_parent)){ + if(!is_alpha_blocker || *blocker_parent!=alpha_parent){ + Simplex_handle sigma_parent(*blocker_parent); + + sigma_parent.difference(alpha_parent); + + if (sigma_parent.dimension()>=2 && sigma_parent.first_vertex() == x_parent){ + Root_simplex_handle sigma_id(parent_complex.get_id(sigma_parent)); + auto sigma_link = this->get_simplex_address(sigma_id); + // ie if the vertices of sigma are vertices of the link + if(sigma_link){ + bool is_new_blocker = true; + for(auto a : alpha_parent){ + for(auto eta_parent : parent_complex.const_blocker_range(a)){ + if(!is_alpha_blocker || *eta_parent!=alpha_parent){ + Simplex_handle eta_minus_alpha(*eta_parent); + eta_minus_alpha.difference(alpha_parent); + if(eta_minus_alpha != sigma_parent && + sigma_parent.contains_difference(*eta_parent,alpha_parent)){ + is_new_blocker = false; + break; + } + } + } + if (!is_new_blocker) break; + } + if (is_new_blocker) + this->add_blocker(new Simplex_handle(*sigma_link)); + + } + } + } + } + } + } + + +public: + /** + * @brief compute vertices, edges and blockers of the link. + * @details If the boolean only_superior_vertices is true, then the link is computed only + * with vertices that are greater than vertices of alpha_parent_adress. + */ + void build_link(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool is_alpha_blocker = false) + { + assert(is_alpha_blocker || parent_complex.contains(alpha_parent_adress)); + compute_link_vertices(parent_complex,alpha_parent_adress,only_superior_vertices_); + compute_link_edges(parent_complex,alpha_parent_adress,is_alpha_blocker); + compute_link_blockers(parent_complex,alpha_parent_adress,is_alpha_blocker); + } + + + /** + * @brief build the link of a blocker which is the link + * of the blocker's simplex if this simplex had been + * removed from the blockers of the complex. + */ + friend void build_link_of_blocker( + const ComplexType & parent_complex, + Simplex_handle& blocker, + Skeleton_blocker_link_complex & result){ + assert(blocker.dimension()>=2); + assert(parent_complex.contains_blocker(blocker)); + result.clear(); + result.build_link(parent_complex,blocker,true); + } + + +}; + +} + +} // namespace GUDHI + +#endif /* SKELETON_BLOCKERS_LINK_COMPLEX_H_ */ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h new file mode 100644 index 00000000..4f9060ef --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h @@ -0,0 +1,525 @@ +/* + * Skeleton_blocker_simplifiable_complex.h + * + * Created on: Feb 4, 2014 + * Author: dsalinas + */ + +#ifndef GUDHI_SKELETON_BLOCKERS_SIMPLIFIABLE_COMPLEX_H_ +#define GUDHI_SKELETON_BLOCKERS_SIMPLIFIABLE_COMPLEX_H_ + +#include "gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h" + +namespace Gudhi{ + +namespace skbl { + +/** + * \brief Class that allows simplification operation on a simplicial complex represented + * by a skeleton/blockers pair. + */ +template<typename SkeletonBlockerDS> +class Skeleton_blocker_simplifiable_complex : public Skeleton_blocker_complex<SkeletonBlockerDS> +{ + template<class ComplexType> friend class Skeleton_blocker_sub_complex; + +public: + + typedef Skeleton_blocker_complex<SkeletonBlockerDS> SkeletonBlockerComplex; + + typedef typename SkeletonBlockerComplex::Graph_edge Graph_edge; + + typedef typename SkeletonBlockerComplex::boost_adjacency_iterator boost_adjacency_iterator; + typedef typename SkeletonBlockerComplex::Edge_handle Edge_handle; + typedef typename SkeletonBlockerComplex::boost_vertex_handle boost_vertex_handle; + typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle; + typedef typename SkeletonBlockerComplex::Root_vertex_handle Root_vertex_handle; + typedef typename SkeletonBlockerComplex::Simplex_handle Simplex_handle; + typedef typename SkeletonBlockerComplex::Root_simplex_handle Root_simplex_handle; + typedef typename SkeletonBlockerComplex::Blocker_handle Blocker_handle; + + + typedef typename SkeletonBlockerComplex::Root_simplex_iterator Root_simplex_iterator; + typedef typename SkeletonBlockerComplex::Simplex_handle_iterator Simplex_handle_iterator; + typedef typename SkeletonBlockerComplex::BlockerMap BlockerMap; + typedef typename SkeletonBlockerComplex::BlockerPair BlockerPair; + typedef typename SkeletonBlockerComplex::BlockerMapIterator BlockerMapIterator; + typedef typename SkeletonBlockerComplex::BlockerMapConstIterator BlockerMapConstIterator; + + typedef typename SkeletonBlockerComplex::Visitor Visitor; + + + /** @name Constructors / Destructors / Initialization + */ + //@{ + Skeleton_blocker_simplifiable_complex(int num_vertices_ = 0,Visitor* visitor_=NULL): + Skeleton_blocker_complex<SkeletonBlockerDS>(num_vertices_,visitor_){ } + + + /** + * @brief Constructor with a list of simplices + * @details The list of simplices must be the list + * of simplices of a simplicial complex, sorted with increasing dimension. + * todo take iterator instead + */ + Skeleton_blocker_simplifiable_complex(std::list<Simplex_handle>& simplices,Visitor* visitor_=NULL): + Skeleton_blocker_complex<SkeletonBlockerDS>(simplices,visitor_) + {} + + + + virtual ~Skeleton_blocker_simplifiable_complex(){ + } + + //@} + + /** + * Returns true iff the blocker 'sigma' is popable. + * To define popable, let us call 'L' the complex that + * consists in the current complex without the blocker 'sigma'. + * A blocker 'sigma' is then "popable" if the link of 'sigma' + * in L is reducible. + * + */ + virtual bool is_popable_blocker(Blocker_handle sigma) const{ + assert(this->contains_blocker(*sigma)); + Skeleton_blocker_link_complex<Skeleton_blocker_simplifiable_complex> link_blocker_sigma; + build_link_of_blocker(*this,*sigma,link_blocker_sigma); + + bool res = link_blocker_sigma.is_contractible()==CONTRACTIBLE; + return res; + } + + + + +private: + /** + * @returns the list of blockers of the simplex + * + * @todo a enlever et faire un iterateur sur tous les blockers a la place + */ + std::list<Blocker_handle> get_blockers(){ + std::list<Blocker_handle> res; + for (auto blocker : this->blocker_range()){ + res.push_back(blocker); + } + return res; + } + + + +public: + + /** + * Removes all the popable blockers of the complex and delete them. + * @returns the number of popable blockers deleted + */ + void remove_popable_blockers(){ + std::list<Vertex_handle> vertex_to_check; + for(auto v : this->vertex_range()) + vertex_to_check.push_front(v); + + while(!vertex_to_check.empty()){ + Vertex_handle v = vertex_to_check.front(); + vertex_to_check.pop_front(); + + bool blocker_popable_found=true; + while (blocker_popable_found){ + blocker_popable_found = false; + for(auto block : this->blocker_range(v)){ + if (this->is_popable_blocker(block)) { + for(Vertex_handle nv : *block) + if(nv!=v) vertex_to_check.push_back(nv); + this->delete_blocker(block); + blocker_popable_found = true; + break; + } + } + } + } + } + + + /** + * Removes all the popable blockers of the complex passing through v and delete them. + */ + void remove_popable_blockers(Vertex_handle v){ + bool blocker_popable_found=true; + while (blocker_popable_found){ + blocker_popable_found = false; + for(auto block : this->blocker_range(v)){ + if (is_popable_blocker(block)) { + this->delete_blocker(block); + blocker_popable_found = true; + } + } + } + } + + + /** + * Remove the star of the vertex 'v' + */ + void remove_star(Vertex_handle v){ + // we remove the blockers that are not consistent anymore + + update_blockers_after_remove_star_of_vertex_or_edge(v); + + while (this->degree(v) > 0) + { + Vertex_handle w( * (adjacent_vertices(v.vertex, this->skeleton).first)); + this->remove_edge(v,w); + } + this->remove_vertex(v); + } + +private: + /** + * after removing the star of a simplex, blockers sigma that contains this simplex must be removed. + * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added + * whenever the dimension of tau is at least 2. + */ + void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed){ + std::list <Blocker_handle> blockers_to_update; + if(simplex_to_be_removed.empty()) return; + + auto v0 = simplex_to_be_removed.first_vertex(); + for (auto blocker : this->blocker_range(v0)){ + if(blocker->contains(simplex_to_be_removed)) + blockers_to_update.push_back(blocker); + } + + for(auto blocker_to_update : blockers_to_update){ + Simplex_handle sub_blocker_to_be_added; + bool sub_blocker_need_to_be_added = + (blocker_to_update->dimension()-simplex_to_be_removed.dimension()) >= 2; + if(sub_blocker_need_to_be_added){ + sub_blocker_to_be_added = *blocker_to_update; + sub_blocker_to_be_added.difference(simplex_to_be_removed); + } + this->delete_blocker(blocker_to_update); + if(sub_blocker_need_to_be_added) + this->add_blocker(sub_blocker_to_be_added); + } + } + + + +public: + /** + * Remove the star of the edge connecting vertices a and b. + * @returns the number of blocker that have been removed + */ + void remove_star(Vertex_handle a, Vertex_handle b){ + update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a,b)); + // we remove the edge + this->remove_edge(a,b); + } + + + /** + * Remove the star of the edge 'e'. + */ + void remove_star(Edge_handle e){ + return remove_star(this->first_vertex(e),this->second_vertex(e)); + } + + /** + * Remove the star of the simplex 'sigma' which needs to belong to the complex + */ + void remove_star(const Simplex_handle& sigma){ + assert(this->contains(sigma)); + if (sigma.dimension()==0) + remove_star(sigma.first_vertex()); + else + if (sigma.dimension()==1) + remove_star(sigma.first_vertex(),sigma.last_vertex()); + else + update_blockers_after_remove_star_of_simplex(sigma); + } + +private: + void update_blockers_after_remove_star_of_simplex(const Simplex_handle& sigma){ + std::list <Blocker_handle> blockers_to_remove; + for (auto blocker : this->blocker_range(sigma.first_vertex())){ + if(blocker->contains(sigma)) + blockers_to_remove.push_back(blocker); + } + for(auto blocker_to_update : blockers_to_remove) + this->delete_blocker(blocker_to_update); + this->add_blocker(sigma); + } + +public: + + enum simplifiable_status{ NOT_HOMOTOPY_EQ,MAYBE_HOMOTOPY_EQ,HOMOTOPY_EQ}; + simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex){ + return MAYBE_HOMOTOPY_EQ; + } + + + + enum contractible_status{ NOT_CONTRACTIBLE,MAYBE_CONTRACTIBLE,CONTRACTIBLE}; + /** + * @brief %Test if the complex is reducible using a strategy defined in the class + * (by default it tests if the complex is a cone) + * @details Note that NO could be returned if some invariant ensures that the complex + * is not a point (for instance if the euler characteristic is different from 1). + * This function will surely have to return MAYBE in some case because the + * associated problem is undecidable but it in practice, it can often + * be solved with the help of geometry. + */ + virtual contractible_status is_contractible() const{ + if (this->is_cone()) return CONTRACTIBLE; + else return MAYBE_CONTRACTIBLE; + // return this->is_cone(); + } + + + /** @Edge contraction operations + */ + //@{ + + + + /** + * @return If ignore_popable_blockers is true + * then the result is true iff the link condition at edge ab is satisfied + * or equivalently iff no blocker contains ab. + * If ignore_popable_blockers is false then the + * result is true iff all blocker containing ab are popable. + */ + bool link_condition(Vertex_handle a, Vertex_handle b,bool ignore_popable_blockers = false) const{ + for (auto blocker : this->const_blocker_range(a)) + if ( blocker->contains(b) ){ + // false if ignore_popable_blockers is false + // otherwise the blocker has to be popable + return ignore_popable_blockers && is_popable_blocker(blocker); + } + return true; + } + + /** + * @return If ignore_popable_blockers is true + * then the result is true iff the link condition at edge ab is satisfied + * or equivalently iff no blocker contains ab. + * If ignore_popable_blockers is false then the + * result is true iff all blocker containing ab are popable. + */ + bool link_condition(Edge_handle & e,bool ignore_popable_blockers = false) const{ + const Graph_edge& edge = (*this)[e]; + assert(this->get_address(edge.first())); + assert(this->get_address(edge.second())); + Vertex_handle a(*this->get_address(edge.first())); + Vertex_handle b(*this->get_address(edge.second())); + return link_condition(a,b,ignore_popable_blockers); + } + + +protected: + /** + * Compute simplices beta such that a.beta is an order 0 blocker + * that may be used to construct a new blocker after contracting ab. + * Suppose that the link condition Link(ab) = Link(a) inter Link(b) + * is satisfied. + */ + void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const{ + for (auto const & blocker : this->const_blocker_range(a)) + { + Simplex_handle beta = (*blocker); + beta.remove_vertex(a); + buffer.push_back(beta); + } + + Simplex_handle n; + this->add_neighbours(b,n); + this->remove_neighbours(a,n); + n.remove_vertex(a); + + + for (Vertex_handle y : n) + { + Simplex_handle beta; + beta.add_vertex( y ); + buffer.push_back(beta); + } + } + + +private: + + /** + * @brief "Replace" the edge 'bx' by the edge 'ax'. + * Assume that the edge 'bx' was present whereas 'ax' was not. + * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'. + * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added + * and just before edge 'bx' had been removed. That way, it can + * eventually access to information of 'ax'. + */ + void swap_edge(Vertex_handle a,Vertex_handle b,Vertex_handle x){ + this->add_edge(a,x); + if (this->visitor) this->visitor->on_swaped_edge(a,b,x); + this->remove_edge(b,x); + } + + +private: + /** + * @brief removes all blockers passing through the edge 'ab' + */ + void delete_blockers_around_vertex(Vertex_handle v){ + std::list <Blocker_handle> blockers_to_delete; + for (auto blocker : this->blocker_range(v)){ + blockers_to_delete.push_back(blocker); + } + while (!blockers_to_delete.empty()){ + this->remove_blocker(blockers_to_delete.back()); + blockers_to_delete.pop_back(); + } + + } + /** + * @brief removes all blockers passing through the edge 'ab' + */ + void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b){ + std::list<Blocker_handle> blocker_to_delete; + for (auto blocker : this->blocker_range(a)) + if (blocker->contains(b)) blocker_to_delete.push_back(blocker); + while (!blocker_to_delete.empty()) + { + this->delete_blocker(blocker_to_delete.back()); + blocker_to_delete.pop_back(); + } + } + +public: + + /** + * Contracts the edge. + * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, + * it removes first all blockers passing through 'ab' + */ + void contract_edge(Edge_handle edge){ + contract_edge(this->first_vertex(edge),this->second_vertex(edge)); + } + + + /** + * Contracts the edge connecting vertices a and b. + * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, + * it removes first all blockers passing through 'ab' + */ + void contract_edge(Vertex_handle a, Vertex_handle b){ + assert(this->contains_vertex(a)); + assert(this->contains_vertex(b)); + assert(this->contains_edge(a,b)); + + // if some blockers passes through 'ab', we remove them. + if (!link_condition(a,b)) + delete_blockers_around_edge(a,b); + + std::set<Simplex_handle> blockers_to_add; + + get_blockers_to_be_added_after_contraction(a,b,blockers_to_add); + + delete_blockers_around_vertices(a,b); + + update_edges_after_contraction(a,b); + + this->remove_vertex(b); + + notify_changed_edges(a); + + for(auto block : blockers_to_add) + this->add_blocker(block); + + assert(this->contains_vertex(a)); + assert(!this->contains_vertex(b)); + } + + +private: + + void get_blockers_to_be_added_after_contraction(Vertex_handle a,Vertex_handle b,std::set<Simplex_handle>& blockers_to_add){ + blockers_to_add.clear(); + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType; + + LinkComplexType link_a(*this,a); + LinkComplexType link_b(*this,b); + + std::vector<Simplex_handle> vector_alpha, vector_beta; + + tip_blockers(a,b,vector_alpha); + tip_blockers(b,a,vector_beta); + + for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha){ + for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) + { + Simplex_handle sigma = *alpha; sigma.union_vertices(*beta); + Root_simplex_handle sigma_id = this->get_id(sigma); + if ( this->contains(sigma) && + proper_faces_in_union<SkeletonBlockerComplex>(sigma_id,link_a,link_b)) + { + // Blocker_handle blocker = new Simplex_handle(sigma); + sigma.add_vertex(a); + blockers_to_add.insert(sigma); + } + } + } + } + + + /** + * delete all blockers that passes through a or b + */ + void delete_blockers_around_vertices(Vertex_handle a,Vertex_handle b){ + std::vector<Blocker_handle> blocker_to_delete; + for(auto bl : this->blocker_range(a)) + blocker_to_delete.push_back(bl); + for(auto bl : this->blocker_range(b)) + blocker_to_delete.push_back(bl); + while (!blocker_to_delete.empty()) + { + this->delete_blocker(blocker_to_delete.back()); + blocker_to_delete.pop_back(); + } + } + + + void update_edges_after_contraction(Vertex_handle a,Vertex_handle b){ + // We update the set of edges + this->remove_edge(a,b); + + // For all edges {b,x} incident to b, + // we remove {b,x} and add {a,x} if not already there. + while (this->degree(b)> 0) + { + Vertex_handle x(*(adjacent_vertices(b.vertex, this->skeleton).first)); + if(!this->contains_edge(a,x)) + // we 'replace' the edge 'bx' by the edge 'ax' + this->swap_edge(a,b,x); + else + this->remove_edge(b,x); + } + } + + void notify_changed_edges(Vertex_handle a){ + // We notify the visitor that all edges incident to 'a' had changed + boost_adjacency_iterator v, v_end; + + for (tie(v, v_end) = adjacent_vertices(a.vertex, this->skeleton); v != v_end; ++v) + if (this->visitor) this->visitor->on_changed_edge(a,Vertex_handle(*v)); + } + + //@} + + +}; + +} + +} // namespace GUDHI + +#endif /* GUDHI_SKELETON_BLOCKERS_SIMPLIFIABLE_COMPLEX_H_ */ diff --git a/src/Skeleton_blocker/test/CMakeLists.txt b/src/Skeleton_blocker/test/CMakeLists.txt new file mode 100644 index 00000000..d66dcf64 --- /dev/null +++ b/src/Skeleton_blocker/test/CMakeLists.txt @@ -0,0 +1,7 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIskbl) + +add_executable ( TestSkeletonBlockerComplex TestSkeletonBlockerComplex.cpp ) +add_executable ( TestSimplifiable TestSimplifiable.cpp ) + + diff --git a/src/Skeleton_blocker/test/TestSimplifiable.cpp b/src/Skeleton_blocker/test/TestSimplifiable.cpp new file mode 100644 index 00000000..9f802463 --- /dev/null +++ b/src/Skeleton_blocker/test/TestSimplifiable.cpp @@ -0,0 +1,281 @@ +/* + * TestSimplifiable.cxx + * + * Created on: Feb 4, 2014 + * Author: dsalinas + */ + + +#include <stdio.h> +#include <stdlib.h> +#include <string> +#include <fstream> +#include <sstream> +#include "gudhi/Test.h" +//#include "Skeleton_blocker/Simplex.h" +#include "gudhi/Skeleton_blocker_complex.h" +#include "gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h" +#include "gudhi/Skeleton_blocker_simplifiable_complex.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h" + + +using namespace std; + +using namespace Gudhi; + +using namespace skbl; + +template<typename ComplexType> class Skeleton_blocker_sub_complex; +typedef Skeleton_blocker_simplifiable_complex<Skeleton_blocker_simple_traits> Complex; +typedef Complex::Vertex_handle Vertex_handle; +typedef Complex::Root_vertex_handle Root_vertex_handle; +typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle; +// true iff v \in complex +bool assert_vertex(Complex &complex,Vertex_handle v){ + Simplex_handle simplex(v); + bool test = complex.contains(simplex); + assert(test); + return test; +} + +// true iff the blocker (a,b,c) is in complex +bool assert_blocker(Complex &complex,Root_vertex_handle a,Root_vertex_handle b,Root_vertex_handle c){ + + return complex.contains_blocker(Simplex_handle(*complex.get_address(a),*complex.get_address(b),*complex.get_address(c))); + //return complex.contains_blocker((a),(b),(c)); +} + +void build_complete(int n,Complex& complex){ + complex.clear(); + for(int i=0;i<n;i++) + complex.add_vertex(); + for(int i=0;i<n;i++) + for(int j=0;j<i;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); +} + +bool test_contraction1(){ + enum { a, b, x, y, z, n }; + Complex complex(n); + build_complete(n,complex); + complex.remove_edge(Vertex_handle(b),Vertex_handle(z)); + complex.add_blocker(Vertex_handle(a),Vertex_handle(x),Vertex_handle(y)); + complex.add_blocker(Vertex_handle(b),Vertex_handle(x),Vertex_handle(y)); + + // Print result + cerr << "complex before complex"<< complex.to_string()<<endl; + + cerr <<endl<<endl; + complex.contract_edge(Vertex_handle(a),Vertex_handle(b)); + // Print result + cerr << "ContractEdge(0,1)\n"; + PRINT(complex.to_string()); + + // verification + for (int i=0;i<5;i++) + if (i!=1) assert_vertex(complex,Vertex_handle(i)); + bool test1 = !complex.contains_edge(Vertex_handle(a),Vertex_handle(b)); + bool test2 = assert_blocker(complex,Root_vertex_handle(a),Root_vertex_handle(x),Root_vertex_handle(y)); + bool test3 = complex.num_edges()==6; + bool test4 = complex.num_blockers()==1; + Simplex_handle sigma; + sigma.add_vertex(Vertex_handle(a)); + sigma.add_vertex(Vertex_handle(x)); + sigma.add_vertex(Vertex_handle(y)); + sigma.add_vertex(Vertex_handle(z)); + bool test5 = !(complex.contains(sigma)); + return test1&&test2&&test3&&test4&&test5; +} + + +bool test_contraction2(){ + enum { a, b, x, y, z, n }; + Complex complex(n); + build_complete(n,complex); + complex.remove_edge(Vertex_handle(b),Vertex_handle(x)); + Simplex_handle blocker; + blocker.add_vertex(Vertex_handle(a)); + blocker.add_vertex(Vertex_handle(y)); + blocker.add_vertex(Vertex_handle(z)); + + complex.add_blocker(blocker); + + // Print result + cerr << "complex complex"<< complex.to_string(); + cerr <<endl<<endl; + complex.contract_edge(Vertex_handle(a),Vertex_handle(b)); + + cerr << "complex.ContractEdge(a,b)"<< complex.to_string(); + + cerr <<endl<<endl; + + // there should be one blocker (a,c,d,e) in the complex + bool test ; + test = complex.contains_blocker(Simplex_handle(a,x,y,z)); + test = test && complex.num_blockers()==1; + return test; +} + +bool test_link_condition1(){ + + Complex complex(0); + // Build the complexes + build_complete(4,complex); + complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2)); + + + // Print result + cerr << "complex complex"<< complex.to_string(); + cerr <<endl<<endl; + + bool weak_link_condition = complex.link_condition(Vertex_handle(1),Vertex_handle(2),true); + + bool strong_link_condition = complex.link_condition(Vertex_handle(1),Vertex_handle(2),false); + + return weak_link_condition && !strong_link_condition; +} + + + + +bool test_collapse1(){ + // xxx implement remove_star(simplex) before + + Complex complex(5); + build_complete(4,complex); + complex.add_vertex(); + complex.add_edge(2,4); + complex.add_edge(3,4); + // Print result + cerr << "initial complex :\n"<< complex.to_string(); + cerr <<endl<<endl; + + Simplex_handle simplex_123(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3)); + complex.remove_star(simplex_123); + cerr << "complex.remove_star(1,2,3):\n"<< complex.to_string(); + cerr <<endl<<endl; + + // verification + bool blocker123_here = complex.contains_blocker(simplex_123); + cerr <<"----> Ocomplex \n"; + return blocker123_here; +} + +bool test_collapse2(){ + Complex complex(5); + build_complete(4,complex); + complex.add_vertex(); + complex.add_edge(Vertex_handle(1),Vertex_handle(4)); + complex.add_edge(Vertex_handle(2),Vertex_handle(4)); + complex.add_edge(Vertex_handle(3),Vertex_handle(4)); + complex.add_blocker(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3),Vertex_handle(4)); + // Print result + cerr << "initial complex :\n"<< complex.to_string(); + cerr <<endl<<endl; + + Simplex_handle sigma(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3)); + complex.remove_star(sigma); + cerr << "complex.remove_star(1,2,3):\n"<< complex.to_string(); + cerr <<endl<<endl; + + // verification + bool blocker_removed = !complex.contains_blocker(Simplex_handle(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3),Vertex_handle(4))); + bool blocker123_here = complex.contains_blocker(sigma); + return blocker_removed && blocker123_here; +} + +bool test_collapse3(){ + Complex complex(5); + build_complete(4,complex); + complex.add_vertex(); + complex.add_edge(Vertex_handle(1),Vertex_handle(4)); + complex.add_edge(Vertex_handle(2),Vertex_handle(4)); + complex.add_edge(Vertex_handle(3),Vertex_handle(4)); + complex.add_blocker(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3),Vertex_handle(4)); + // Print result + cerr << "initial complex:\n"<< complex.to_string(); + cerr <<endl<<endl; + + complex.remove_star(Vertex_handle(2)); + cerr << "complex after remove star of 2:\n"<< complex.to_string(); + + bool blocker134_here = complex.contains_blocker(Simplex_handle(Vertex_handle(1),Vertex_handle(3),Vertex_handle(4))); + bool blocker1234_here = complex.contains_blocker(Simplex_handle(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3),Vertex_handle(4))); + // verification + // assert_blocker(complex,1,2,3); + // assert(!complex.ContainsBlocker(new AddressSimplex(1,2,3,4))); + cerr <<"----> Ocomplex \n"; + return blocker134_here && !blocker1234_here; +} + + + +bool test_remove_popable_blockers(){ + Complex complex; + build_complete(4,complex); + complex.add_vertex(); + complex.add_edge(Vertex_handle(3),Vertex_handle(4)); + complex.add_edge(Vertex_handle(2),Vertex_handle(4)); + Simplex_handle sigma1=Simplex_handle(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3)); + Simplex_handle sigma2=Simplex_handle(Vertex_handle(2),Vertex_handle(3),Vertex_handle(4)); + + complex.add_blocker(sigma1); + complex.add_blocker(sigma2); + cerr << "complex complex"<< complex.to_string(); + cerr <<endl<<endl; + cerr << "complex.RemovePopableBlockers();"<<endl; + complex.remove_popable_blockers(); + cerr << "complex complex"<< complex.to_string(); + cerr <<endl<<endl; + + bool test1 = (complex.num_blockers()==1); + + + // test 2 + complex.clear(); + build_complete(4,complex); + complex.add_vertex(); + complex.add_vertex(); + complex.add_edge(Vertex_handle(3),Vertex_handle(4)); + complex.add_edge(Vertex_handle(2),Vertex_handle(4)); + complex.add_edge(Vertex_handle(2),Vertex_handle(5)); + complex.add_edge(Vertex_handle(3),Vertex_handle(5)); + complex.add_edge(Vertex_handle(4),Vertex_handle(5)); + sigma1=Simplex_handle(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3)); + sigma2=Simplex_handle(Vertex_handle(2),Vertex_handle(3),Vertex_handle(4)); + + complex.add_blocker(sigma1); + complex.add_blocker(sigma2); + cerr << "complex complex"<< complex.to_string(); + cerr <<endl<<endl; cerr << "complex.RemovePopableBlockers();"<<endl; + complex.remove_popable_blockers(); + cerr << "complex complex"<< complex.to_string(); + + cerr <<endl<<endl; + bool test2 = (complex.num_blockers()==0); + return test1&&test2; +} + + + +int main (int argc, char *argv[]) +{ + Tests tests_simplifiable_complex; + tests_simplifiable_complex.add("Test contraction 1",test_contraction1); + tests_simplifiable_complex.add("Test contraction 2",test_contraction2); + tests_simplifiable_complex.add("Test Link condition 1",test_link_condition1); + tests_simplifiable_complex.add("Test remove popable blockers",test_remove_popable_blockers); + + + tests_simplifiable_complex.add("Test collapse 1",test_collapse1); + tests_simplifiable_complex.add("Test collapse 2",test_collapse2); + tests_simplifiable_complex.add("Test collapse 3",test_collapse3); + + + tests_simplifiable_complex.run(); + + if(tests_simplifiable_complex.run()) + return EXIT_SUCCESS; + else + return EXIT_FAILURE; +} diff --git a/src/Skeleton_blocker/test/TestSkeletonBlockerComplex.cpp b/src/Skeleton_blocker/test/TestSkeletonBlockerComplex.cpp new file mode 100644 index 00000000..c89f84e2 --- /dev/null +++ b/src/Skeleton_blocker/test/TestSkeletonBlockerComplex.cpp @@ -0,0 +1,773 @@ +#include <stdio.h> +#include <stdlib.h> +#include <string> +#include <fstream> +#include <sstream> +#include "gudhi/Utils.h" +#include "gudhi/Test.h" +//#include "Skeleton_blocker/Simplex.h" +#include "gudhi/Skeleton_blocker_complex.h" +#include "gudhi/Skeleton_blocker_link_complex.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h" +//#include "Simple_vertex.h" +//#include "Simple_edge.h" + +using namespace std; + +using namespace Gudhi; + +using namespace skbl; + + +typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex; +typedef Complex::Vertex_handle Vertex_handle; +typedef Complex::Root_vertex_handle Root_vertex_handle; +typedef Complex::Simplex_handle Simplex_handle; +typedef Complex::Root_simplex_handle Root_simplex_handle; +typedef Simplex_handle::Simplex_vertex_const_iterator Simplex_vertex_const_iterator; +typedef Complex::Edge_handle Edge_handle; + +// true iff v \in complex +bool assert_vertex(Complex &complex,Vertex_handle v){ + //assert(complex.contains(v)); + return complex.contains(v); +} + +bool assert_simplex(Complex &complex,Root_vertex_handle a,Root_vertex_handle b,Root_vertex_handle c){ + return true; + // AddressSimplex simplex((a),(b),(c)); + // return complex.contains(&simplex); +} + +// true iff the blocker (a,b,c) is in complex +bool assert_blocker(Complex &complex,Root_vertex_handle a,Root_vertex_handle b,Root_vertex_handle c){ + return true; + //return complex.contains_blocker((a),(b),(c)); +} + +// true iff the blocker (a,b,c,d) is in complex +bool assert_blocker(Complex &complex,Root_vertex_handle a,Root_vertex_handle b,Root_vertex_handle c,Root_vertex_handle d){ + return true; + //Simplex blocker (a,b,c,d); + //return complex.contains_blocker(&blocker); +} + + +void build_complete(int n,Complex& complex){ + complex.clear(); + for(int i=0;i<n;i++) + complex.add_vertex(); + +// for(int i=n-1;i>=0;i--) +// for(int j=i-1;j>=0;j--) +// complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + + for(int i=0;i<n;i++) + for(int j=0;j<i;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); +} + + +bool test_simplex(){ +// PRINT("test simplex"); + Simplex_handle simplex(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2),Vertex_handle(3)); + for (auto i = simplex.begin() ; i != simplex.end() ; ++i){ + PRINT(*i); + auto j = i; + for (++j ; + j != simplex.end() ; + ++j){ + PRINT(*j); + } + } + return simplex.dimension()==3; +} + + +bool test_iterator_vertices1(){ + int n = 10; + Complex complex(10); + cerr << "complex.num_vertices():"<<complex.num_vertices()<<endl; + int num_vertex_seen = 0; + for(auto vi :complex.vertex_range()){ + cerr << "vertex:"<<vi<<endl; + ++num_vertex_seen; + } + return num_vertex_seen == n; +} + +bool test_iterator_vertices2(){ + int n = 10; + Complex complex(10); + build_complete(10,complex); + cerr << "complex.num_vertices():"<<complex.num_vertices()<<endl; + cerr << "complex.num_edges():"<<complex.num_edges()<<endl; + int num_vertex_seen = 0; + for(auto vi :complex.vertex_range(Vertex_handle(2))){ + cerr << "vertex:"<<vi<<endl; + ++num_vertex_seen; + } + std::cerr<<"num_vertex_seen:"<<num_vertex_seen<<std::endl; + return num_vertex_seen == (n-1); +} + + + +bool test_iterator_edge(){ + const int n = 10; + Complex complex(n); + for(int i=0;i<n;i++) + for(int j=0;j<i;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + complex.remove_edge(Vertex_handle(2),Vertex_handle(3)); + complex.remove_edge(Vertex_handle(3),Vertex_handle(5)); + cerr << "complex.num_edges():"<<complex.num_edges()<<endl; + int num_edges_seen = 0; + for(auto edge : complex.edge_range()){ + cerr << "edge :"<<complex[edge]<<endl; + ++num_edges_seen; + } + + return num_edges_seen == n*(n-1)/2-2; +} + +bool test_iterator_edge2(){ + const int n = 10; + Complex complex(n); + for(int i=0;i<n;i++) + for(int j=0;j<i;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + complex.remove_edge(Vertex_handle(2),Vertex_handle(3)); + complex.remove_edge(Vertex_handle(3),Vertex_handle(5)); + cerr << "complex.num_edges():"<<complex.num_edges()<<endl; + int num_neigbors_seen = 0; + for(auto neighbor : complex.vertex_range(Vertex_handle(2))){ + cerr << "neighbor"<<neighbor<<endl; + ++num_neigbors_seen; + } + return num_neigbors_seen==8; +} + + + +bool test_iterator_edge3(){ + const int n = 10; + Complex complex(n); + for(int i=0;i<n;i++) + for(int j=0;j<i;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + complex.remove_edge(Vertex_handle(2),Vertex_handle(3)); + complex.remove_edge(Vertex_handle(3),Vertex_handle(5)); + cerr << "complex.num_edges():"<<complex.num_edges()<<endl; + int num_neigbors_seen = 0; + for(auto edge : complex.edge_range(Vertex_handle(2))){ + std::cerr << edge<< std::endl; + ++num_neigbors_seen; + } + return num_neigbors_seen==8; +} + + + +bool test_iterator_triangles(){ + const int n = 7; + Complex complex(n); + //create a "ring" around '0' + for(int i=1;i<n;i++) + complex.add_edge(Vertex_handle(0),Vertex_handle(i)); + for(int i=1;i<n-1;i++) + complex.add_edge(Vertex_handle(i),Vertex_handle(i+1)); + complex.add_edge(Vertex_handle(1),Vertex_handle(6)); + + PRINT(complex.to_string()); + + int num_triangles_seen=0; + //for (auto t : complex.triangle_range(5)){ + TEST("triangles around 5 (should be 2 of them):"); + for (auto t : complex.triangle_range(5)){ + PRINT(t); + ++num_triangles_seen; + } + bool test = (num_triangles_seen==2); + + num_triangles_seen=0; + TEST("triangles around 0 (should be 6 of them):"); + for (auto t : complex.triangle_range(0)){ + PRINT(t); + ++num_triangles_seen; + } + test = test&&(num_triangles_seen==6); + + // we now add another triangle + complex.add_vertex(); + complex.add_edge(Vertex_handle(4),Vertex_handle(7)); + complex.add_edge(Vertex_handle(3),Vertex_handle(7)); + complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(6)); + num_triangles_seen=0; + + TEST("triangles (should be 6 of them):"); + num_triangles_seen=0; + for (auto t : complex.triangle_range()){ + PRINT(t); + ++num_triangles_seen; + } + test = test&&(num_triangles_seen==6); + PRINT(num_triangles_seen); + + return test; +} + + +//#include "combinatorics/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h" + +bool test_iterator_simplices(){ + Complex complex(6); + complex.add_edge(Vertex_handle(0),Vertex_handle(1)); + complex.add_edge(Vertex_handle(1),Vertex_handle(2)); + complex.add_edge(Vertex_handle(2),Vertex_handle(0)); + complex.add_edge(Vertex_handle(1),Vertex_handle(3)); + complex.add_edge(Vertex_handle(2),Vertex_handle(3)); + complex.add_edge(Vertex_handle(2),Vertex_handle(5)); + complex.add_edge(Vertex_handle(3),Vertex_handle(5)); + complex.add_edge(Vertex_handle(2),Vertex_handle(4)); + complex.add_edge(Vertex_handle(4),Vertex_handle(5)); + complex.add_edge(Vertex_handle(3),Vertex_handle(4)); + + complex.add_blocker(Vertex_handle(2),Vertex_handle(3),Vertex_handle(4),Vertex_handle(5)); + + bool correct_number_simplices = true; + + std::map<Vertex_handle,unsigned> expected_num_simplices; + + expected_num_simplices[Vertex_handle(0)] = 4; + expected_num_simplices[Vertex_handle(1)] = 6; + expected_num_simplices[Vertex_handle(2)] = 11; + expected_num_simplices[Vertex_handle(3)] = 9; + expected_num_simplices[Vertex_handle(4)] = 7; + expected_num_simplices[Vertex_handle(5)] = 7; + + for(auto pair : expected_num_simplices){ + unsigned num_simplices_around = 0; + for(const auto& simplex : complex.simplex_range(pair.first)){ + simplex.dimension(); + DBGVALUE(simplex); + ++num_simplices_around; + } + + correct_number_simplices = correct_number_simplices && (num_simplices_around == pair.second); + + DBGMSG("current vertex:",pair.first); + DBGMSG("expected_num_simplices:",pair.second); + DBGMSG("found:",num_simplices_around); + } + return correct_number_simplices; +} + + + +bool test_iterator_simplices2(){ + Complex complex(2); + complex.add_edge(Vertex_handle(0),Vertex_handle(1)); + + for(const auto& s:complex.triangle_range()){ + s.dimension(); + return false; // there are no triangles + } + + unsigned num_simplices = 0 ; + + + DBGVALUE(complex.to_string()); + + for(const auto& simplex : complex.simplex_range(Vertex_handle(0))){ + simplex.dimension(); + DBGVALUE(simplex); + } + + + for(const auto& simplex : complex.simplex_range()){ + DBGVALUE(simplex); + simplex.dimension(); + ++num_simplices; + } + bool correct_number_simplices = (num_simplices == 3); + return correct_number_simplices; +} + + +bool test_iterator_simplices3(){ + Complex complex(3); + complex.add_edge(Vertex_handle(0),Vertex_handle(1)); + complex.add_edge(Vertex_handle(1),Vertex_handle(2)); + complex.add_edge(Vertex_handle(2),Vertex_handle(0)); + complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2)); + + unsigned num_simplices = 0 ; + + for(const auto& simplex : complex.simplex_range(Vertex_handle(0))){ + simplex.dimension(); + DBGVALUE(simplex); + } + + + for(const auto& simplex : complex.simplex_range()){ + DBGVALUE(simplex); + simplex.dimension(); + ++num_simplices; + } + bool correct_number_simplices = (num_simplices == 6); + return correct_number_simplices; +} + +bool test_iterator_simplices4(){ + Complex empty_complex; + for(auto v : empty_complex.vertex_range()){ + v; + } + for(auto e : empty_complex.edge_range()){ + empty_complex[e]; + } + for(auto t : empty_complex.triangle_range()){ + t.dimension(); + } + for(auto s : empty_complex.simplex_range()){ + s.dimension(); + } + return true; +} + + + + + +template<typename Map> +auto blocker_range(Map map) -> decltype( map | boost::adaptors::map_values){ + return map| boost::adaptors::map_values ; +} + + +bool test_iterator_blockers(){ + Complex complex; + Simplex_handle alpha; + Simplex_handle vertex_set_expected; + // Build the complexes + for (int i=0;i<20;i++){ + complex.add_vertex(); + } + for (int i=10;i<15;i++){ + for (int j=i+1;j<15;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + } + + complex.add_blocker(Simplex_handle(Vertex_handle(10),Vertex_handle(11),Vertex_handle(12))); + complex.add_blocker(Simplex_handle(Vertex_handle(2),Vertex_handle(1),Vertex_handle(10))); + complex.add_blocker(Simplex_handle(Vertex_handle(10),Vertex_handle(9),Vertex_handle(15))); + complex.add_blocker(Simplex_handle(Vertex_handle(1),Vertex_handle(9),Vertex_handle(8))); + + // Print result + int num_blockers=0; + for(auto blockers : complex.blocker_range(Vertex_handle(10))){ + TESTVALUE(*blockers) ; + num_blockers++; + } + bool test = (num_blockers==3); + + num_blockers=0; + for (auto blockers : complex.blocker_range()){ + TESTVALUE(*blockers) ; + num_blockers++; + } + test = test && (num_blockers==4) ; + + return test; +} + + +bool test_link0(){ + + enum { a, b, c, d, n }; + Complex complex(n); + complex.add_edge(Vertex_handle(b),Vertex_handle(c));complex.add_edge(Vertex_handle(c),Vertex_handle(d)); + Simplex_handle alpha = Simplex_handle(Vertex_handle(c)); + Skeleton_blocker_link_complex<Complex> L(complex,alpha); + PRINT(L.num_vertices()); + PRINT(L.to_string()); + + bool test1 = L.contains_vertex(*L.get_address(Root_vertex_handle(b))); + bool test2 = L.contains_vertex(*L.get_address(Root_vertex_handle(d))); + bool test3 = L.num_edges()==0; + bool test4 = L.num_blockers()==0; + return test1&&test2&&test3&&test4; + +} + +bool test_link1(){ + Complex complex; + + + // Build the complexes + for (int i=0;i<20;i++){ + complex.add_vertex(); + } + for (int i=10;i<15;i++){ + for (int j=i+1;j<15;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + } + Simplex_handle alpha(Vertex_handle(12),Vertex_handle(14)); + Skeleton_blocker_link_complex<Complex> L(complex,alpha); + // Complexes built + + // verification + bool test1 = L.contains_vertex(*L.get_address(Root_vertex_handle(10))); + bool test2 = L.contains_vertex(*L.get_address(Root_vertex_handle(11))); + bool test3 = L.contains_vertex(*L.get_address(Root_vertex_handle(13))); + bool test4 = L.num_edges()==3; + bool test5 = L.num_blockers()==0; + Root_simplex_handle simplex; + simplex.add_vertex(Root_vertex_handle(10)); + simplex.add_vertex(Root_vertex_handle(11)); + simplex.add_vertex(Root_vertex_handle(13)); + bool test6(L.get_simplex_address(simplex)); + bool test7 = L.contains(*(L.get_simplex_address(simplex))); + cerr <<"----> Ocomplex \n"; + return test1&&test2&&test3&&test4&&test5&&test6&&test7 ; + +} + + +bool test_link2(){ + Complex complex; + + Simplex_handle alpha; + Simplex_handle vertex_set_expected; + // Build the complexes + for (int i=0;i<20;i++){ + complex.add_vertex(); + } + for (int i=10;i<15;i++){ + for (int j=i+1;j<15;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + } + complex.add_blocker(Vertex_handle(10),Vertex_handle(11),Vertex_handle(13)); + alpha = Simplex_handle(Vertex_handle(12),Vertex_handle(14)); + Skeleton_blocker_link_complex<Complex> L(complex,alpha); + // Complexes built + + // Print result + cerr << "complex complex"<< complex.to_string(); + cerr <<endl<<endl; + cerr << "L= Link_complex("<<alpha<<") : \n"<<L.to_string(); + + // verification + bool test1 = L.contains_vertex(*L.get_address(Root_vertex_handle(10))); + bool test2 = L.contains_vertex(*L.get_address(Root_vertex_handle(11))); + bool test3 = L.contains_vertex(*L.get_address(Root_vertex_handle(13))); + bool test4 = L.num_edges()==3; + bool test5 = L.num_blockers()==1; + Root_simplex_handle simplex; + simplex.add_vertex(Root_vertex_handle(10)); + simplex.add_vertex(Root_vertex_handle(11)); + simplex.add_vertex(Root_vertex_handle(13)); + bool test6 = L.contains_blocker(*(L.get_simplex_address(simplex))); + cerr <<"----> Ocomplex \n"; + return test1&&test2&&test3&&test4&&test5&&test6 ; +} + +bool test_link3(){ + Complex complex; + + Simplex_handle alpha; + Simplex_handle vertex_set_expected; + // Build the complexes + for (int i=0;i<20;i++){ + complex.add_vertex(); + } + for (int i=10;i<15;i++){ + for (int j=i+1;j<15;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + } + complex.add_blocker(Vertex_handle(10),Vertex_handle(11),Vertex_handle(12)); + alpha = Simplex_handle(Vertex_handle(12),Vertex_handle(14)); + Skeleton_blocker_link_complex<Complex> L(complex,alpha); + // Complexes built + + // Print result + cerr << "complex complex"<< complex.to_string(); + cerr <<endl<<endl; + cerr << "L= Link_complex("<<alpha<<") : \n"<<L.to_string(); + + // verification + bool test = assert_vertex(L,*L.get_address(Root_vertex_handle(10))); + test = test&& assert_vertex(L,*L.get_address(Root_vertex_handle(11))); + test = test&& assert_vertex(L,*L.get_address(Root_vertex_handle(13))); + test = test&& L.num_edges()==2; + test = test&&L.contains_edge(*L.get_address(Root_vertex_handle(10)),*L.get_address(Root_vertex_handle(13))); + test=test&&L.contains_edge(*L.get_address(Root_vertex_handle(13)),*L.get_address(Root_vertex_handle(11))); + test=test&&L.num_blockers()==0; + return test; +} + +bool test_link4(){ + Complex complex; + + // Build the complexes + for (int i=0;i<20;i++){ + complex.add_vertex(); + } + for (int i=10;i<15;i++){ + for (int j=i+1;j<15;j++) + complex.add_edge(Vertex_handle(i),Vertex_handle(j)); + } + complex.add_blocker(Vertex_handle(10),Vertex_handle(11),Vertex_handle(12),Vertex_handle(13)); + Simplex_handle alpha(Vertex_handle(12),Vertex_handle(14)); + Skeleton_blocker_link_complex<Complex> L(complex,alpha); + // Complexes built + + // verification + bool test1 = L.contains_vertex(*L.get_address(Root_vertex_handle(10))); + bool test2 = L.contains_vertex(*L.get_address(Root_vertex_handle(11))); + bool test3 = L.contains_vertex(*L.get_address(Root_vertex_handle(13))); + bool test4 = L.num_edges()==3; + bool test5 = L.num_blockers()==1; + Root_simplex_handle simplex; + simplex.add_vertex(Root_vertex_handle(10)); + simplex.add_vertex(Root_vertex_handle(11)); + simplex.add_vertex(Root_vertex_handle(13)); + bool test6 = L.contains_blocker(*(L.get_simplex_address(simplex))); + cerr <<"----> Ocomplex \n"; + return test1&&test2&&test3&&test4&&test5&&test6 ; + +} + +bool test_link5(){ + Complex complex(0,new Print_complex_visitor<Vertex_handle>()); + // Build the complexes + build_complete(4,complex); + complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2),Vertex_handle(3)); + + Simplex_handle alpha(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2)); + + + Skeleton_blocker_link_complex<Complex> L(complex,alpha); // Complexes built + + // Print result + PRINT(complex.to_string()); + cerr <<endl<<endl; + PRINT(L.to_string()); + + // verification + return L.num_vertices()==0; +} + +bool test_link6(){ + Complex complex(0,new Print_complex_visitor<Vertex_handle>()); + // Build the complexes + build_complete(4,complex); + complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2)); + + Simplex_handle alpha(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2)); + + Skeleton_blocker_link_complex<Complex> link_blocker_alpha; + + build_link_of_blocker(complex,alpha,link_blocker_alpha); + + // Print result + PRINT(complex.to_string()); + cerr <<endl<<endl; + PRINT(link_blocker_alpha.to_string()); + + // verification + return link_blocker_alpha.num_vertices()==1; +} + + +bool test_link7(){ + Complex complex(0,new Print_complex_visitor<Vertex_handle>()); + // Build the complexes + build_complete(6,complex); + complex.add_vertex(); + complex.add_vertex(); + for(int i = 3; i<6; ++i){ + complex.add_edge(Vertex_handle(i),Vertex_handle(6)); + complex.add_edge(Vertex_handle(i),Vertex_handle(7)); + } + complex.add_edge(Vertex_handle(6),Vertex_handle(7)); + complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2)); + complex.add_blocker(Vertex_handle(3),Vertex_handle(4),Vertex_handle(5)); + + Simplex_handle alpha(Vertex_handle(3),Vertex_handle(4),Vertex_handle(5)); + + Skeleton_blocker_link_complex<Complex> link_blocker_alpha; + + build_link_of_blocker(complex,alpha,link_blocker_alpha); + + //the result should be the edge {6,7} plus the blocker {0,1,2} + + // Print result + PRINT(complex.to_string()); + cerr <<endl<<endl; + DBGVALUE(link_blocker_alpha.to_string()); + + Skeleton_blocker_link_complex<Complex> link_blocker_alpha_cpy = link_blocker_alpha; + + DBGVALUE(link_blocker_alpha_cpy.to_string()); + + bool equal_complexes = + (link_blocker_alpha.num_vertices() == link_blocker_alpha_cpy.num_vertices()) + &&(link_blocker_alpha.num_blockers() == link_blocker_alpha_cpy.num_blockers()) + &&(link_blocker_alpha.num_edges() == link_blocker_alpha_cpy.num_edges()) + ; + DBGVALUE((link_blocker_alpha.num_blockers() == link_blocker_alpha_cpy.num_blockers())); + DBGVALUE((link_blocker_alpha.num_blockers() )); + DBGVALUE(( link_blocker_alpha_cpy.num_blockers())); + + DBGVALUE(equal_complexes); + + // verification + return link_blocker_alpha.num_vertices()==5 && link_blocker_alpha.num_edges()==4 && link_blocker_alpha.num_blockers()==1 && equal_complexes; +} + + + + + + + + +template<typename SimplexHandle> +void add_triangle_edges(int a,int b,int c,list<SimplexHandle>& simplices){ + typedef SimplexHandle Simplex_handle; + typedef typename SimplexHandle::Vertex_handle Vertex_handle; + + simplices.push_back(Simplex_handle(Vertex_handle(a),Vertex_handle(b) )); + simplices.push_back(Simplex_handle(Vertex_handle(b),Vertex_handle(c) )); + simplices.push_back(Simplex_handle(Vertex_handle(c),Vertex_handle(a) )); +} + +template<typename SimplexHandle> +void add_triangle(int a,int b,int c,list<SimplexHandle>& simplices){ + typedef SimplexHandle Simplex_handle; + typedef typename SimplexHandle::Vertex_handle Vertex_handle; + simplices.push_back(Simplex_handle(Vertex_handle(a),Vertex_handle(b),Vertex_handle(c))); +} + +bool test_constructor(){ + list <Simplex_handle> simplices; + + simplices.push_back(Simplex_handle(Vertex_handle(0))); + simplices.push_back(Simplex_handle(Vertex_handle(1))); + simplices.push_back(Simplex_handle(Vertex_handle(2))); + simplices.push_back(Simplex_handle(Vertex_handle(3))); + simplices.push_back(Simplex_handle(Vertex_handle(4))); + simplices.push_back(Simplex_handle(Vertex_handle(5))); + + simplices.push_back(Simplex_handle(Vertex_handle(3),Vertex_handle(5) )); + + add_triangle_edges(0,1,5,simplices); + add_triangle_edges(1,2,3,simplices); + add_triangle_edges(2,3,4,simplices); + add_triangle_edges(1,3,4,simplices); + add_triangle_edges(1,2,4,simplices); + + + add_triangle(0,1,5,simplices); + add_triangle(1,2,3,simplices); + add_triangle(2,3,4,simplices); + add_triangle(1,3,4,simplices); + add_triangle(1,2,4,simplices); + + + Complex complex(simplices); + + PRINT(complex.to_string()); + + return ( complex.num_vertices()==6&&complex.num_edges()==10&& complex.num_blockers()==2); +} + + +list<Simplex_handle> subfaces(Simplex_handle top_face){ + list<Simplex_handle> res; + if(top_face.dimension()==-1) return res; + if(top_face.dimension()==0) { + res.push_back(top_face); + return res; + } + else{ + Vertex_handle first_vertex = top_face.first_vertex(); + top_face.remove_vertex(first_vertex); + res = subfaces(top_face); + list<Simplex_handle> copy = res; + for(auto& simplex : copy){ + simplex.add_vertex(first_vertex); + } + res.push_back(Simplex_handle(first_vertex)); + res.splice(res.end(),copy); + return res; + } +} + + +bool test_constructor2(){ + Simplex_handle simplex; + for(int i =0 ; i < 5;++i) + simplex.add_vertex(i); + + list <Simplex_handle> simplices(subfaces(simplex)); + simplices.remove(simplex); + + Complex complex(simplices); + + PRINT(complex.to_string()); + + for(auto b : complex.const_blocker_range()){ + cout << "b:"<<b<<endl; + } + + return ( complex.num_vertices()==5&&complex.num_edges()==10&& complex.num_blockers()==1); +} + + + + +int main (int argc, char *argv[]) +{ + Tests tests_complex; + tests_complex.add("test simplex",test_simplex); + tests_complex.add("test_link0",test_link0); + tests_complex.add("test_link1",test_link1); + tests_complex.add("test_link2",test_link2); + tests_complex.add("test_link3",test_link3); + tests_complex.add("test_link4",test_link4); + tests_complex.add("test_link5",test_link5); + tests_complex.add("test_link6",test_link6); + tests_complex.add("test_link7",test_link7); + + tests_complex.add("test iterator vertices 1",test_iterator_vertices1); + tests_complex.add("test iterator vertices 2",test_iterator_vertices2); + tests_complex.add("test iterator edges",test_iterator_edge); + tests_complex.add("test iterator edges 2",test_iterator_edge2); + tests_complex.add("test iterator edges 3",test_iterator_edge3); + + tests_complex.add("test iterator simplices",test_iterator_simplices); + tests_complex.add("test iterator simplices2",test_iterator_simplices2); + tests_complex.add("test iterator simplices3",test_iterator_simplices3); + tests_complex.add("test iterator simplices4",test_iterator_simplices4); + + + tests_complex.add("test iterator blockers",test_iterator_blockers); + tests_complex.add("test_iterator_triangles",test_iterator_triangles); + + tests_complex.add("test_constructor_list_simplices",test_constructor); + tests_complex.add("test_constructor_list_simplices2",test_constructor2); + + if(tests_complex.run()){ + return EXIT_SUCCESS; + } + else{ + return EXIT_FAILURE; + } + + // test_iterator_simplices(); +} + |