diff options
Diffstat (limited to 'src')
31 files changed, 2973 insertions, 2198 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b910d477..c209a020 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,17 +2,19 @@ cmake_minimum_required(VERSION 2.6) project(GUDHI) -find_package(Boost REQUIRED COMPONENTS system filesystem unit_test_framework chrono timer REQUIRED) +find_package(Boost REQUIRED COMPONENTS system filesystem program_options chrono timer REQUIRED) if(MSVC) # Turn off some VC++ warnings SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4267 /wd4668 /wd4311 /wd4800 /wd4820 /wd4503 /wd4244 /wd4345 /wd4996 /wd4396 /wd4018") else() - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 -std=c++11 -Wall -Wpedantic -lboost_system -lgmpxx -lgmp") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wpedantic") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g -O0") set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") endif() +# BOOST ISSUE result_of vs C++11 +add_definitions(-DBOOST_RESULT_OF_USE_DECLTYPE) find_package(Boost) if(NOT Boost_FOUND) @@ -22,6 +24,8 @@ else() LINK_DIRECTORIES(${Boost_LIBRARY_DIRS}) include_directories(include/) + add_subdirectory(example/Simplex_tree) + add_subdirectory(example/Persistent_cohomology) add_subdirectory(example/Skeleton_blocker) add_subdirectory(example/Contraction) diff --git a/src/Contraction/example/CMakeLists.txt b/src/Contraction/example/CMakeLists.txt index 717830da..13eb5537 100644 --- a/src/Contraction/example/CMakeLists.txt +++ b/src/Contraction/example/CMakeLists.txt @@ -1,6 +1,7 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIskbl) -add_executable ( RipsContraction Rips_contraction.cpp) -target_link_libraries( RipsContraction ${Boost_TIMER_LIBRARY}) +add_executable(RipsContraction Rips_contraction.cpp) +target_link_libraries(RipsContraction ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) +add_test(RipsContraction ${CMAKE_CURRENT_BINARY_DIR}/RipsContraction ${CMAKE_SOURCE_DIR}/data/meshes/sphere3D_2646.off 0.2) diff --git a/src/Contraction/test/CMakeLists.txt b/src/Contraction/test/CMakeLists.txt new file mode 100644 index 00000000..6dddff1b --- /dev/null +++ b/src/Contraction/test/CMakeLists.txt @@ -0,0 +1,20 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIskbl) + +if(NOT MSVC) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +endif() + +add_executable(TestContraction TestContraction.cpp) + +# Unitary tests +add_test(TestContraction ${CMAKE_CURRENT_BINARY_DIR}/TestContraction ${CMAKE_SOURCE_DIR}/data/meshes/SO3_50000.off 0.2) + +if (LCOV_PATH) + # Lcov code coverage of unitary test + add_test(src/Contraction/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Contraction) +endif() + +cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Contraction/include/gudhi") diff --git a/src/Contraction/test/TestContraction.cpp b/src/Contraction/test/TestContraction.cpp new file mode 100644 index 00000000..872f24a0 --- /dev/null +++ b/src/Contraction/test/TestContraction.cpp @@ -0,0 +1,198 @@ + /* 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-Mediterranee (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/>. + */ +#include <ctime> +#include <list> + +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h" +#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h" +//#include "Skeleton_blocker/Simplex.h" +#include "gudhi/Skeleton_blocker_contractor.h" +#include "gudhi/Utils.h" +#include "gudhi/Test.h" +#include "gudhi/Skeleton_blocker_geometric_complex.h" +//#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> + + +#include "gudhi/Contraction/Edge_profile.h" +#include "gudhi/Contraction/policies/Cost_policy.h" +#include "gudhi/Contraction/policies/Edge_length_cost.h" +#include "gudhi/Contraction/policies/Placement_policy.h" +#include "gudhi/Contraction/policies/Middle_placement.h" +#include "gudhi/Contraction/policies/Valid_contraction_policy.h" +#include "gudhi/Contraction/policies/Dummy_valid_contraction.h" +#include "gudhi/Contraction/policies/Link_condition_valid_contraction.h" + + +using namespace std; + +using namespace Gudhi; + +using namespace skbl; + +struct Geometry_trait{ + typedef std::vector<double> Point; +}; + +typedef Geometry_trait::Point Point; + + +typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> AbstractComplex; +typedef Skeleton_blocker_simple_geometric_traits<Geometry_trait> Complex_geometric_traits; + + +typedef Skeleton_blocker_geometric_complex< Complex_geometric_traits > Complex; + +typedef Complex::Vertex_handle Vertex_handle; +typedef Complex::Simplex_handle Simplex_handle; + +typedef Complex::Root_vertex_handle Root_vertex_handle; + +using namespace contraction; + +typedef Skeleton_blocker_contractor<Complex> Complex_contractor; + +typedef Edge_profile<Complex> Profile; + +// compute the distance todo utiliser Euclidean_geometry a la place +template<typename Point> +double eucl_distance(const Point& a,const Point& b){ + double res = 0; + auto a_coord = a.begin(); + auto b_coord = b.begin(); + for(; a_coord != a.end(); a_coord++, b_coord++){ + res += (*a_coord - *b_coord) * (*a_coord - *b_coord); + } + return sqrt(res); +} + +// build the Rips complex todo utiliser Euclidean_geometry a la place +template<typename ComplexType> +void build_rips(ComplexType& complex, double offset){ + if (offset<=0) return; + auto vertices = complex.vertex_range(); + for (auto p = vertices.begin(); p != vertices.end(); ++p) + for (auto q = p; ++q != vertices.end(); /**/) + if (eucl_distance(complex.point(*p), complex.point(*q)) < 2 * offset){ + complex.add_edge(*p,*q); + } +} + + + + +void test_contraction_rips(string name_file, double offset){ + Complex complex; + // load the points + if (!read_off_file<Complex>(name_file,complex,true)){ + std::cerr << "Unable to read file:"<<name_file<<std::endl; + std::cerr << "current path : "; + system("pwd"); + std::cerr<<endl; + return; + } + + clock_t time = clock(); + + TEST("build the Rips complex"); + + build_rips(complex,offset); + + std::cerr << "Rips contruction took "<< ( (float)(clock()-time))/CLOCKS_PER_SEC << " seconds\n"; + + TESTMSG("Initial number of vertices :",complex.num_vertices()); + TESTMSG("Initial number of edges :",complex.num_edges()); + TESTMSG("Initial number of blockers:",complex.num_blockers()); + + time = clock(); + + Complex_contractor contractor(complex, + new Edge_length_cost<Profile>, + contraction::make_first_vertex_placement<Profile>(), + contraction::make_link_valid_contraction<Profile>(), + contraction::make_remove_popable_blockers_visitor<Profile>()); + contractor.contract_edges(); + + TESTVALUE(complex.to_string()); + + TESTVALUE(complex.num_vertices()); + TESTVALUE(complex.num_edges()); + TESTVALUE(complex.num_blockers()); + + std::cerr << "Edge contractions took "<< ( (float)(clock()-time))/CLOCKS_PER_SEC << " seconds\n"; + +} + + +void test_geometric_link(){ + + Complex complex; + std::vector<double> p0(2,0); + std::vector<double> p1(2,0); p1[0] = 1.; + std::vector<double> p2(2,1); + complex.add_vertex(p0); + complex.add_vertex(p1); + complex.add_vertex(p2); + + 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)); + + + + cerr << "complex points:" <<endl; + for(auto v : complex.vertex_range()){ + cerr <<v <<" -> "; + DBGCONT(complex.point(v)); + } + + cerr << "complex : "<<complex.to_string()<<endl; + + + + + auto link = complex.link(Vertex_handle(0)); + + + cerr << "link of 0 points:" <<endl; + for(auto v : link.vertex_range()){ + cerr <<v <<" -> "; + DBGCONT(link.point(v)); + } + + cerr << "link : "<<link.to_string()<<endl; +} + + + + +int main (int argc, char *argv[]) +{ + if (argc!=3){ + std::cerr << "Usage "<<argv[0]<<" GUDHIPATH/src/data/sphere3D.off 0.1 to load the file GUDHIPATH/src/data/sphere3D.off and contract the Rips complex built with paremeter 0.2.\n"; + return -1; + } + + std::string name_file(argv[1]); + test_contraction_rips(name_file,atof(argv[2])); +} + + diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index dca4a04d..5c337f0b 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -1,22 +1,25 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIExPersCohom) -# NEED TO FIND BOOST NEEDED COMPONENTS TO LINK WITH -find_package(Boost 1.45.0 COMPONENTS system program_options) -message("Boost_lib = ${Boost_LIBRARIES}") +add_executable(rips_persistence rips_persistence.cpp) +target_link_libraries(rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_SYSTEM_LIBRARY}) +add_test(rips_persistence_2 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.25 -d 3 -p 2 -m 100) +add_test(rips_persistence_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.25 -d 3 -p 3 -m 100) -add_executable ( rips_persistence rips_persistence.cpp ) -target_link_libraries(rips_persistence ${Boost_LIBRARIES}) -add_executable ( rips_persistence_from_alpha_shape_3 rips_persistence_from_alpha_shape_3.cpp ) -target_link_libraries(rips_persistence_from_alpha_shape_3 ${Boost_LIBRARIES}) +add_executable(rips_persistence_from_alpha_shape_3 rips_persistence_from_alpha_shape_3.cpp ) +target_link_libraries(rips_persistence_from_alpha_shape_3 ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_SYSTEM_LIBRARY}) if(GMPXX_FOUND AND GMP_FOUND) message("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") message("GMP_LIBRARIES = ${GMP_LIBRARIES}") - add_executable ( rips_multifield_persistence rips_multifield_persistence.cpp ) - target_link_libraries(rips_multifield_persistence ${Boost_LIBRARIES} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + add_executable(rips_multifield_persistence rips_multifield_persistence.cpp ) + target_link_libraries(rips_multifield_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + add_test(rips_multifield_persistence_2_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.25 -d 3 -p 2 -q 3 -m 100) + add_test(rips_multifield_persistence_2_71 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.25 -d 3 -p 2 -q 71 -m 100) add_executable ( performance_rips_persistence performance_rips_persistence.cpp ) - target_link_libraries(performance_rips_persistence ${Boost_LIBRARIES} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + target_link_libraries(performance_rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + add_test(rips_persistence_from_alpha_shape_3_2_0 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_from_alpha_shape_3 ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 2 -m 0) + add_test(rips_persistence_from_alpha_shape_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_from_alpha_shape_3 ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 3 -m 100) endif() diff --git a/src/Persistent_cohomology/example/performance_rips_persistence.cpp b/src/Persistent_cohomology/example/performance_rips_persistence.cpp index 37d41ee1..077c2b07 100644 --- a/src/Persistent_cohomology/example/performance_rips_persistence.cpp +++ b/src/Persistent_cohomology/example/performance_rips_persistence.cpp @@ -31,6 +31,7 @@ #include <chrono> using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; /* Compute the persistent homology of the complex cpx with coefficients in Z/pZ. */ template< typename FilteredComplex> diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp index 95ef8809..2505897e 100644 --- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -30,6 +30,7 @@ #include <boost/program_options.hpp> using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; typedef int Vertex_handle; typedef double Filtration_value; @@ -69,6 +70,8 @@ int main (int argc, char * argv[]) st.insert_graph(prox_graph); // insert the proximity graph in the simplex tree st.expansion( dim_max ); // expand the graph until dimension dim_max + std::cout << "st=" << st << std::endl; + // Sort the simplices in the order of the filtration st.initialize_filtration(); diff --git a/src/Persistent_cohomology/example/rips_persistence.cpp b/src/Persistent_cohomology/example/rips_persistence.cpp index d7927223..f7f9527f 100644 --- a/src/Persistent_cohomology/example/rips_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_persistence.cpp @@ -29,6 +29,7 @@ #include <boost/program_options.hpp> using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; typedef int Vertex_handle; typedef double Filtration_value; diff --git a/src/Persistent_cohomology/example/rips_persistence_from_alpha_shape_3.cpp b/src/Persistent_cohomology/example/rips_persistence_from_alpha_shape_3.cpp index ddaafea2..e886aea7 100644 --- a/src/Persistent_cohomology/example/rips_persistence_from_alpha_shape_3.cpp +++ b/src/Persistent_cohomology/example/rips_persistence_from_alpha_shape_3.cpp @@ -29,6 +29,7 @@ #include <boost/program_options.hpp> using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; typedef int Vertex_handle; typedef double Filtration_value; diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index 2b3efc28..27d64b8a 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -1,37 +1,46 @@ - /* 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): Clément Maria - * - * 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/>. - */ +/* 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): Clément Maria + * + * 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 SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_H_ +#define SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_H_ -#ifndef _PERSISTENCECOMPUTATION_SIMPLEXTREE_ -#define _PERSISTENCECOMPUTATION_SIMPLEXTREE_ +#include <gudhi/Persistent_cohomology/Persistent_cohomology_column.h> +#include <gudhi/Persistent_cohomology/Field_Zp.h> #include <boost/tuple/tuple.hpp> #include <boost/intrusive/set.hpp> #include <boost/pending/disjoint_sets.hpp> #include <boost/intrusive/list.hpp> #include <boost/pool/object_pool.hpp> -#include "gudhi/Persistent_cohomology/Persistent_cohomology_column.h" -#include "gudhi/Persistent_cohomology/Field_Zp.h" -namespace Gudhi{ +#include <map> +#include <utility> +#include <list> +#include <vector> +#include <set> + +namespace Gudhi { + +namespace persistent_cohomology { /** \defgroup persistent_cohomology Persistent Cohomology * @@ -50,112 +59,107 @@ namespace Gudhi{ composed of three elements: topological spaces, their homology groups and an evolution scheme. - <DT>Topological Spaces:</DT> - Topological spaces are represented by simplicial complexes. - Let \f$V = \{1, \cdots ,|V|\}\f$ be a set of <EM>vertices</EM>. - A <EM>simplex</EM> \f$\sigma\f$ is a subset of vertices - \f$\sigma \subseteq V\f$. A <EM>simplicial complex</EM> \f$\mathbf{K}\f$ - on \f$V\f$ is a collection of simplices \f$\{\sigma\}\f$, - \f$\sigma \subseteq V\f$, such that \f$\tau \subseteq \sigma \in \mathbf{K} - \Rightarrow \tau \in \mathbf{K}\f$. The dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$ - is its number of elements minus 1. A <EM>filtration</EM> of a simplicial complex is - a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq - f(\sigma)\f$ whenever \f$\tau \subseteq \sigma\f$. - - We define the concept FilteredComplex which enumerates the requirements for a class - to represent a filtered complex from which persistent homology may be computed. - We use the vocabulary of simplicial complexes, but the concept - is valid for any type of cell complex. The main requirements - are the definition of: - \li type <CODE>Indexing_tag</CODE>, which is a model of the concept - <CODE>IndexingTag</CODE>, - describing the nature of the indexing scheme, - \li type Simplex_handle to manipulate simplices, - \li method <CODE>int dimension(Simplex_handle)</CODE> returning - the dimension of a simplex, - \li type and method <CODE>Boundary_simplex_range - boundary_simplex_range(Simplex_handle)</CODE> that returns - a range giving access to the codimension 1 subsimplices of the - input simplex, as-well-as the coefficients \f$(-1)^i\f$ in the - definition of the operator \f$\partial\f$. The iterators have - value type <CODE>Simplex_handle</CODE>, - \li type and method - <CODE>Filtration_simplex_range filtration_simplex_range ()</CODE> - that returns a range giving - access to all the simplices of the complex read in the order - assigned by the indexing scheme, - \li type and method - <CODE>Filtration_value filtration (Simplex_handle)</CODE> that returns the value of - the filtration on the simplex represented by the handle. - - <DT>Homology:</DT> - For a ring \f$\mathcal{R}\f$, the group of <EM>n-chains</EM>, - denoted \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$, of \f$\mathbf{K}\f$ is the - group of formal sums of - n-simplices with \f$\mathcal{R}\f$ coefficients. The <EM>boundary operator</EM> is a - linear operator - \f$\partial_n: \mathbf{C}_n(\mathbf{K},\mathcal{R}) \rightarrow \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R})\f$ - such that \f$\partial_n \sigma = \partial_n [v_0, \cdots , v_n] = - \sum_{i=0}^n (-1)^{i}[v_0,\cdots ,\widehat{v_i}, \cdots,v_n]\f$, - where \f$\widehat{v_i}\f$ means \f$v_i\f$ is omitted from the list. The chain - groups form a sequence: - - \f[\cdots \ \ \mathbf{C}_n(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_n\ } \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R}) - \xrightarrow{\partial_{n-1}} \cdots \xrightarrow{\ \partial_2 \ } - \mathbf{C}_1(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_1 \ } \mathbf{C}_0(\mathbf{K},\mathcal{R}) \f] - - of finitely many groups \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$ and homomorphisms - \f$\partial_n\f$, indexed by the dimension \f$n \geq 0\f$. - The boundary operators satisfy the property \f$\partial_n \circ \partial_{n+1}=0\f$ - for every \f$n > 0\f$ - and we define the homology groups: - - \f[\mathbf{H}_n(\mathbf{K},\mathcal{R}) = \ker \partial_n / \mathrm{im} \ \partial_{n+1}\f] - - We refer to \cite Munkres-elementsalgtop1984 for an introduction to homology - theory and to \cite DBLP:books/daglib/0025666 for an introduction to persistent homology. - - <DT>Indexing Scheme:</DT> - "Changing" a simplicial complex consists in applying a simplicial map. - An <EM>indexing scheme</EM> is a directed graph together with a traversal - order, such that two - consecutive nodes in the graph are connected by an arrow (either forward or backward). - The nodes represent simplicial complexes and the directed edges simplicial maps. - - From the computational point of view, there are two types of indexing schemes of - interest - in persistent homology: <EM>linear</EM> ones - \f$\bullet \longrightarrow \bullet \longrightarrow \cdots \longrightarrow \bullet - \longrightarrow \bullet\f$ - in persistent homology \cite DBLP:journals/dcg/ZomorodianC05 , - and <EM>zigzag</EM> ones - \f$\bullet \longrightarrow \bullet \longleftarrow \cdots - \longrightarrow \bullet - \longleftarrow \bullet \f$ in zigzag persistent - homology \cite DBLP:journals/focm/CarlssonS10. - These indexing schemes have a natural left-to-right traversal order, and we - describe them with ranges and iterators. - In the current release of the Gudhi library, only the linear case is implemented. - - In the following, we consider the case where the indexing scheme is induced - by a filtration. - Ordering the simplices - by increasing filtration values (breaking ties so as a simplex appears after - its subsimplices of same filtration value) provides an indexing scheme. - - - - <DT>Implementations:</DT> - We use the <EM>Compressed Annotation Matrix</EM> of \cite DBLP:conf/esa/BoissonnatDM13 to - implement the - persistent cohomology algorithm of \cite DBLP:journals/dcg/SilvaMV11 - and \cite DBLP:conf/compgeom/DeyFW14 for persistence in the class Persistent_cohomology. - - The coefficient fields available as models of CoefficientField are Field_Zp - for \f$\mathbb{Z}_p\f$ (for any prime p) and Multi_field for the multi-field persistence algorithm - -- computing persistence simultaneously in various coefficient fields -- described - in \cite boissonnat:hal-00922572. - + The theory of homology consists in attaching to a topological space a sequence of + (homology) groups, + capturing global topological features + like connected components, holes, cavities, etc. Persistent homology studies the evolution + -- birth, life and death -- of + these features when the topological space is changing. Consequently, the theory is essentially + composed of three elements: + topological spaces, their homology groups and an evolution scheme. + + <DT>Topological Spaces:</DT> + Topological spaces are represented by simplicial complexes. + Let \f$V = \{1, \cdots ,|V|\}\f$ be a set of <EM>vertices</EM>. + A <EM>simplex</EM> \f$\sigma\f$ is a subset of vertices + \f$\sigma \subseteq V\f$. A <EM>simplicial complex</EM> \f$\mathbf{K}\f$ + on \f$V\f$ is a collection of simplices \f$\{\sigma\}\f$, + \f$\sigma \subseteq V\f$, such that \f$\tau \subseteq \sigma \in \mathbf{K} + \Rightarrow \tau \in \mathbf{K}\f$. The dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$ + is its number of elements minus 1. A <EM>filtration</EM> of a simplicial complex is + a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq + f(\sigma)\f$ whenever \f$\tau \subseteq \sigma\f$. + + We define the concept FilteredComplex which enumerates the requirements for a class + to represent a filtered complex from which persistent homology may be computed. + We use the vocabulary of simplicial complexes, but the concept + is valid for any type of cell complex. The main requirements + are the definition of: + \li type <CODE>Indexing_tag</CODE>, which is a model of the concept + <CODE>IndexingTag</CODE>, + describing the nature of the indexing scheme, + \li type Simplex_handle to manipulate simplices, + \li method <CODE>int dimension(Simplex_handle)</CODE> returning + the dimension of a simplex, + \li type and method <CODE>Boundary_simplex_range + boundary_simplex_range(Simplex_handle)</CODE> that returns + a range giving access to the codimension 1 subsimplices of the + input simplex, as-well-as the coefficients \f$(-1)^i\f$ in the + definition of the operator \f$\partial\f$. The iterators have + value type <CODE>Simplex_handle</CODE>, + \li type and method + <CODE>Filtration_simplex_range filtration_simplex_range ()</CODE> + that returns a range giving + access to all the simplices of the complex read in the order + assigned by the indexing scheme, + \li type and method + <CODE>Filtration_value filtration (Simplex_handle)</CODE> that returns the value of + the filtration on the simplex represented by the handle. + + <DT>Homology:</DT> + For a ring \f$\mathcal{R}\f$, the group of <EM>n-chains</EM>, + denoted \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$, of \f$\mathbf{K}\f$ is the + group of formal sums of + n-simplices with \f$\mathcal{R}\f$ coefficients. The <EM>boundary operator</EM> is a + linear operator + \f$\partial_n: \mathbf{C}_n(\mathbf{K},\mathcal{R}) \rightarrow \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R})\f$ + such that \f$\partial_n \sigma = \partial_n [v_0, \cdots , v_n] = + \sum_{i=0}^n (-1)^{i}[v_0,\cdots ,\widehat{v_i}, \cdots,v_n]\f$, + where \f$\widehat{v_i}\f$ means \f$v_i\f$ is omitted from the list. The chain + groups form a sequence: + + \f[\cdots \ \ \mathbf{C}_n(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_n\ } \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R}) + \xrightarrow{\partial_{n-1}} \cdots \xrightarrow{\ \partial_2 \ } + \mathbf{C}_1(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_1 \ } \mathbf{C}_0(\mathbf{K},\mathcal{R}) \f] + + of finitely many groups \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$ and homomorphisms + \f$\partial_n\f$, indexed by the dimension \f$n \geq 0\f$. + The boundary operators satisfy the property \f$\partial_n \circ \partial_{n+1}=0\f$ + for every \f$n > 0\f$ + and we define the homology groups: + + \f[\mathbf{H}_n(\mathbf{K},\mathcal{R}) = \ker \partial_n / \mathrm{im} \ \partial_{n+1}\f] + + We refer to \cite Munkres-elementsalgtop1984 for an introduction to homology + theory and to \cite DBLP:books/daglib/0025666 for an introduction to persistent homology. + + <DT>Indexing Scheme:</DT> + "Changing" a simplicial complex consists in applying a simplicial map. + An <EM>indexing scheme</EM> is a directed graph together with a traversal + order, such that two + consecutive nodes in the graph are connected by an arrow (either forward or backward). + The nodes represent simplicial complexes and the directed edges simplicial maps. + + From the computational point of view, there are two types of indexing schemes of + interest + in persistent homology: <EM>linear</EM> ones + \f$\bullet \longrightarrow \bullet \longrightarrow \cdots \longrightarrow \bullet + \longrightarrow \bullet\f$ + in persistent homology \cite DBLP:journals/dcg/ZomorodianC05 , + and <EM>zigzag</EM> ones + \f$\bullet \longrightarrow \bullet \longleftarrow \cdots + \longrightarrow \bullet + \longleftarrow \bullet \f$ in zigzag persistent + homology \cite DBLP:journals/focm/CarlssonS10. + These indexing schemes have a natural left-to-right traversal order, and we + describe them with ranges and iterators. + In the current release of the Gudhi library, only the linear case is implemented. + + In the following, we consider the case where the indexing scheme is induced + by a filtration. + Ordering the simplices + by increasing filtration values (breaking ties so as a simplex appears after + its subsimplices of same filtration value) provides an indexing scheme. \section Examples We provide several example files: run these examples with -h for details on their use, and read the README file. @@ -180,606 +184,579 @@ different representations. */ /** \brief Computes the persistent cohomology of a filtered complex. -* -* The computation is implemented with a Compressed Annotation Matrix -* (CAM)\cite DBLP:conf/esa/BoissonnatDM13, -* and is adapted to the computation of Multi-Field Persistent Homology (MF) -* \cite boissonnat:hal-00922572 . -* -* \implements PersistentHomology -* -*/ -//Memory allocation policy: classic, use a mempool, etc.*/ - template < class FilteredComplex - , class CoefficientField - > //to do mem allocation policy: classic, mempool, etc. - class Persistent_cohomology { + * + * The computation is implemented with a Compressed Annotation Matrix + * (CAM)\cite DBLP:conf/esa/BoissonnatDM13, + * and is adapted to the computation of Multi-Field Persistent Homology (MF) + * \cite boissonnat:hal-00922572 . + * + * \implements PersistentHomology + * + */ +// Memory allocation policy: classic, use a mempool, etc.*/ +template<class FilteredComplex, class CoefficientField> // to do mem allocation policy: classic, mempool, etc. +class Persistent_cohomology { public: - - typedef FilteredComplex Complex_ds; + typedef FilteredComplex Complex_ds; // Data attached to each simplex to interface with a Property Map. - typedef typename Complex_ds::Simplex_key Simplex_key; - typedef typename Complex_ds::Simplex_handle Simplex_handle; - typedef typename Complex_ds::Filtration_value Filtration_value; - typedef typename CoefficientField::Element Arith_element; + typedef typename Complex_ds::Simplex_key Simplex_key; + typedef typename Complex_ds::Simplex_handle Simplex_handle; + typedef typename Complex_ds::Filtration_value Filtration_value; + typedef typename CoefficientField::Element Arith_element; // Compressed Annotation Matrix types: // Column type - typedef Persistent_cohomology_column < Simplex_key - , Arith_element - > Column; // contains 1 set_hook + typedef Persistent_cohomology_column<Simplex_key, Arith_element> Column; // contains 1 set_hook // Cell type - typedef typename Column::Cell Cell; // contains 2 list_hooks + typedef typename Column::Cell Cell; // contains 2 list_hooks // Remark: constant_time_size must be false because base_hook_cam_h has auto_unlink link_mode - typedef boost::intrusive::list < Cell - , boost::intrusive::constant_time_size<false> - , boost::intrusive::base_hook< base_hook_cam_h > - > Hcell; - - typedef boost::intrusive::set < Column - , boost::intrusive::constant_time_size<false> - > Cam; + typedef boost::intrusive::list<Cell, + boost::intrusive::constant_time_size<false>, + boost::intrusive::base_hook<base_hook_cam_h> > Hcell; + + typedef boost::intrusive::set<Column, + boost::intrusive::constant_time_size<false> > Cam; // Sparse column type for the annotation of the boundary of an element. - typedef std::vector< std::pair<Simplex_key - , Arith_element > > A_ds_type; + typedef std::vector<std::pair<Simplex_key, Arith_element> > A_ds_type; // Persistent interval type. The Arith_element field is used for the multi-field framework. - typedef boost::tuple< Simplex_handle - , Simplex_handle - , Arith_element > Persistent_interval; - -/** \brief Initializes the Persistent_cohomology class. - * - * @param[in] cpx Complex for which the persistent homology is compiuted. - cpx is a model of FilteredComplex - * - * @param[in] persistence_dim_max if true, the persistent homology for the maximal dimension in the - * complex is computed. If false, it is ignored. Default is false. - */ -Persistent_cohomology ( Complex_ds & cpx - , bool persistence_dim_max = false ) -: cpx_ (&cpx) -, dim_max_(cpx.dimension()) // upper bound on the dimension of the simplices -, coeff_field_() // initialize the field coefficient structure. -, ds_rank_ (cpx_->num_simplices()) // union-find -, ds_parent_(cpx_->num_simplices()) // union-find -, ds_repr_ (cpx_->num_simplices(),NULL) // union-find -> annotation vectors -, dsets_(&ds_rank_[0],&ds_parent_[0]) // union-find -, cam_() // collection of annotation vectors -, zero_cocycles_() // union-find -> Simplex_key of creator for 0-homology -, transverse_idx_() // key -> row -, persistent_pairs_() -, interval_length_policy(&cpx,0) -, column_pool_(new boost::object_pool< Column > ()) // memory pools for the CAM -, cell_pool_(new boost::object_pool< Cell > ()) -{ - if( persistence_dim_max ) { ++dim_max_; } - Simplex_key idx_fil = 0; - for(auto & sh : cpx_->filtration_simplex_range()) - { - cpx_->assign_key(sh,idx_fil); ++idx_fil; - dsets_.make_set(cpx_->key(sh)); + typedef boost::tuple<Simplex_handle, Simplex_handle, Arith_element> Persistent_interval; + + /** \brief Initializes the Persistent_cohomology class. + * + * @param[in] cpx Complex for which the persistent homology is compiuted. + cpx is a model of FilteredComplex + * + * @param[in] persistence_dim_max if true, the persistent homology for the maximal dimension in the + * complex is computed. If false, it is ignored. Default is false. + */ + explicit Persistent_cohomology(Complex_ds& cpx) + : cpx_(&cpx), + dim_max_(cpx.dimension()), // upper bound on the dimension of the simplices + coeff_field_(), // initialize the field coefficient structure. + ds_rank_(cpx_->num_simplices()), // union-find + ds_parent_(cpx_->num_simplices()), // union-find + ds_repr_(cpx_->num_simplices(), NULL), // union-find -> annotation vectors + dsets_(&ds_rank_[0], &ds_parent_[0]), // union-find + cam_(), // collection of annotation vectors + zero_cocycles_(), // union-find -> Simplex_key of creator for 0-homology + transverse_idx_(), // key -> row + persistent_pairs_(), + interval_length_policy(&cpx, 0), + column_pool_(new boost::object_pool<Column>()), // memory pools for the CAM + cell_pool_(new boost::object_pool<Cell>()) { + Simplex_key idx_fil = 0; + for (auto & sh : cpx_->filtration_simplex_range()) { + cpx_->assign_key(sh, idx_fil); + ++idx_fil; + dsets_.make_set(cpx_->key(sh)); + } } -} - -~Persistent_cohomology() -{ -//Clean the remaining columns in the matrix. - for(auto & cam_ref : cam_) { cam_ref.col_.clear(); } -//Clean the transversal lists - for( auto & transverse_ref : transverse_idx_ ) - { transverse_ref.second.row_->clear(); delete transverse_ref.second.row_; } -//Clear the memory pools - delete column_pool_; - delete cell_pool_; -} - -private: -struct length_interval { - length_interval ( Complex_ds * cpx - , Filtration_value min_length) - : cpx_(cpx) - , min_length_(min_length) {} - - bool operator()(Simplex_handle sh1, Simplex_handle sh2) - { return cpx_->filtration(sh2) - cpx_->filtration(sh1) > min_length_; } - - void set_length(Filtration_value new_length) { min_length_ = new_length; } - - Complex_ds * cpx_; - Filtration_value min_length_; -}; - -public: -/** \brief Initializes the coefficient field.*/ -void init_coefficients( int charac ) { coeff_field_.init(charac); } -/** \brief Initializes the coefficient field for multi-field persistent homology.*/ -void init_coefficients( int charac_min - , int charac_max ) { coeff_field_.init(charac_min,charac_max); } + /** \brief Initializes the Persistent_cohomology class. + * + * @param[in] cpx Complex for which the persistent homology is compiuted. + cpx is a model of FilteredComplex + * + * @param[in] persistence_dim_max if true, the persistent homology for the maximal dimension in the + * complex is computed. If false, it is ignored. Default is false. + */ + Persistent_cohomology(Complex_ds& cpx, bool persistence_dim_max) + : Persistent_cohomology(cpx) { + if (persistence_dim_max) { + ++dim_max_; + } + } -/** \brief Compute the persistent homology of the filtered simplicial - * complex. - * - * @param[in] min_interval_length the computation disgards all intervals of length - * less or equal than min_interval_length - * - * Assumes that the filtration provided by the simplicial complex is - * valid. Undefined behavior otherwise. */ -void compute_persistent_cohomology ( Filtration_value min_interval_length = 0 ) -{ - interval_length_policy.set_length(min_interval_length); - // Compute all finite intervals - for( auto sh : cpx_->filtration_simplex_range() ) - { - int dim_simplex = cpx_->dimension(sh); - switch(dim_simplex) { - case 0 : break; - case 1 : update_cohomology_groups_edge( sh ) ; break; - default: update_cohomology_groups( sh, dim_simplex ); break; + ~Persistent_cohomology() { +// Clean the remaining columns in the matrix. + for (auto & cam_ref : cam_) { + cam_ref.col_.clear(); } +// Clean the transversal lists + for (auto & transverse_ref : transverse_idx_) { + transverse_ref.second.row_->clear(); + delete transverse_ref.second.row_; + } +// Clear the memory pools + delete column_pool_; + delete cell_pool_; } - // Compute infinite intervals of dimension 0 - Simplex_key key; - for(auto v_sh : cpx_->skeleton_simplex_range(0)) //for all 0-dimensional simplices - { - key = cpx_->key(v_sh); - - if( ds_parent_[key] == key //root of its tree - && zero_cocycles_.find(key) == zero_cocycles_.end() ) - { - persistent_pairs_.push_back( Persistent_interval ( cpx_->simplex(key) - , cpx_->null_simplex() - , coeff_field_.characteristic() ) - ); + + private: + struct length_interval { + length_interval(Complex_ds * cpx, Filtration_value min_length) + : cpx_(cpx), + min_length_(min_length) { } + + bool operator()(Simplex_handle sh1, Simplex_handle sh2) { + return cpx_->filtration(sh2) - cpx_->filtration(sh1) > min_length_; + } + + void set_length(Filtration_value new_length) { + min_length_ = new_length; + } + + Complex_ds * cpx_; + Filtration_value min_length_; + }; + + public: + /** \brief Initializes the coefficient field.*/ + void init_coefficients(uint16_t charac) { + coeff_field_.init(charac); } - for( auto zero_idx : zero_cocycles_ ) - { - persistent_pairs_.push_back( Persistent_interval ( cpx_->simplex(zero_idx.second) - , cpx_->null_simplex() - , coeff_field_.characteristic() ) - ); + /** \brief Initializes the coefficient field for multi-field persistent homology.*/ + void init_coefficients(uint16_t charac_min, uint16_t charac_max) { + coeff_field_.init(charac_min, charac_max); } -// Compute infinite interval of dimension > 0 - for(auto cocycle : transverse_idx_) - { - persistent_pairs_.push_back( Persistent_interval ( cpx_->simplex (cocycle.first) - , cpx_->null_simplex() - , cocycle.second.characteristics_ ) ); + + /** \brief Compute the persistent homology of the filtered simplicial + * complex. + * + * @param[in] min_interval_length the computation disgards all intervals of length + * less or equal than min_interval_length + * + * Assumes that the filtration provided by the simplicial complex is + * valid. Undefined behavior otherwise. */ + void compute_persistent_cohomology(Filtration_value min_interval_length = 0) { + interval_length_policy.set_length(min_interval_length); + // Compute all finite intervals + for (auto sh : cpx_->filtration_simplex_range()) { + int dim_simplex = cpx_->dimension(sh); + switch (dim_simplex) { + case 0: + break; + case 1: + update_cohomology_groups_edge(sh); + break; + default: + update_cohomology_groups(sh, dim_simplex); + break; + } + } + // Compute infinite intervals of dimension 0 + Simplex_key key; + for (auto v_sh : cpx_->skeleton_simplex_range(0)) { // for all 0-dimensional simplices + key = cpx_->key(v_sh); + + if (ds_parent_[key] == key // root of its tree + && zero_cocycles_.find(key) == zero_cocycles_.end()) { + persistent_pairs_.push_back( + Persistent_interval(cpx_->simplex(key), cpx_->null_simplex(), coeff_field_.characteristic())); + } + } + for (auto zero_idx : zero_cocycles_) { + persistent_pairs_.push_back( + Persistent_interval(cpx_->simplex(zero_idx.second), cpx_->null_simplex(), coeff_field_.characteristic())); + } +// Compute infinite interval of dimension > 0 + for (auto cocycle : transverse_idx_) { + persistent_pairs_.push_back( + Persistent_interval(cpx_->simplex(cocycle.first), cpx_->null_simplex(), cocycle.second.characteristics_)); + } } -} - - - -private: -/** \brief Update the cohomology groups under the insertion of an edge. - * - * The 0-homology is maintained with a simple Union-Find data structure, which - * explains the existance of a specific function of edge insertions. */ -void update_cohomology_groups_edge ( Simplex_handle sigma ) -{ - Simplex_handle u,v; - boost::tie(u,v) = cpx_->endpoints(sigma); - - Simplex_key ku = dsets_.find_set( cpx_->key(u) ); - Simplex_key kv = dsets_.find_set( cpx_->key(v) ); - - if(ku != kv ) { // Destroy a connected component - dsets_.link(ku,kv); - // Keys of the simplices which created the connected components containing - // respectively u and v. - Simplex_key idx_coc_u, idx_coc_v; - auto map_it_u = zero_cocycles_.find(ku); - // If the index of the cocycle representing the class is already ku. - if (map_it_u == zero_cocycles_.end()) { idx_coc_u = ku; } - else { idx_coc_u = map_it_u->second; } - - auto map_it_v = zero_cocycles_.find(kv); - // If the index of the cocycle representing the class is already kv. - if (map_it_v == zero_cocycles_.end()) { idx_coc_v = kv; } - else { idx_coc_v = map_it_v->second; } - - if(cpx_->filtration(cpx_->simplex(idx_coc_u)) - < cpx_->filtration(cpx_->simplex(idx_coc_v)) ) // Kill cocycle [idx_coc_v], which is younger. - { - if(interval_length_policy(cpx_->simplex(idx_coc_v),sigma)) { - persistent_pairs_.push_back ( Persistent_interval ( cpx_->simplex(idx_coc_v) - , sigma - , coeff_field_.characteristic() - ) - ); + + private: + /** \brief Update the cohomology groups under the insertion of an edge. + * + * The 0-homology is maintained with a simple Union-Find data structure, which + * explains the existance of a specific function of edge insertions. */ + void update_cohomology_groups_edge(Simplex_handle sigma) { + Simplex_handle u, v; + boost::tie(u, v) = cpx_->endpoints(sigma); + + Simplex_key ku = dsets_.find_set(cpx_->key(u)); + Simplex_key kv = dsets_.find_set(cpx_->key(v)); + + if (ku != kv) { // Destroy a connected component + dsets_.link(ku, kv); + // Keys of the simplices which created the connected components containing + // respectively u and v. + Simplex_key idx_coc_u, idx_coc_v; + auto map_it_u = zero_cocycles_.find(ku); + // If the index of the cocycle representing the class is already ku. + if (map_it_u == zero_cocycles_.end()) { + idx_coc_u = ku; + } else { + idx_coc_u = map_it_u->second; + } + + auto map_it_v = zero_cocycles_.find(kv); + // If the index of the cocycle representing the class is already kv. + if (map_it_v == zero_cocycles_.end()) { + idx_coc_v = kv; + } else { + idx_coc_v = map_it_v->second; } - // Maintain the index of the 0-cocycle alive. - if( kv != idx_coc_v ) { zero_cocycles_.erase( map_it_v ); } - if( kv == dsets_.find_set(kv) ) { - if( ku != idx_coc_u ) { zero_cocycles_.erase( map_it_u ); } + + if (cpx_->filtration(cpx_->simplex(idx_coc_u)) + < cpx_->filtration(cpx_->simplex(idx_coc_v))) { // Kill cocycle [idx_coc_v], which is younger. + if (interval_length_policy(cpx_->simplex(idx_coc_v), sigma)) { + persistent_pairs_.push_back( + Persistent_interval(cpx_->simplex(idx_coc_v), sigma, coeff_field_.characteristic())); + } + // Maintain the index of the 0-cocycle alive. + if (kv != idx_coc_v) { + zero_cocycles_.erase(map_it_v); + } + if (kv == dsets_.find_set(kv)) { + if (ku != idx_coc_u) { + zero_cocycles_.erase(map_it_u); + } zero_cocycles_[kv] = idx_coc_u; } - } - else // Kill cocycle [idx_coc_u], which is younger. - { - if(interval_length_policy(cpx_->simplex(idx_coc_u),sigma)) { - persistent_pairs_.push_back ( Persistent_interval ( cpx_->simplex(idx_coc_u) - , sigma - , coeff_field_.characteristic() - ) - ); + } else { // Kill cocycle [idx_coc_u], which is younger. + if (interval_length_policy(cpx_->simplex(idx_coc_u), sigma)) { + persistent_pairs_.push_back( + Persistent_interval(cpx_->simplex(idx_coc_u), sigma, coeff_field_.characteristic())); } - // Maintain the index of the 0-cocycle alive. - if( ku != idx_coc_u ) { zero_cocycles_.erase( map_it_u ); } - if( ku == dsets_.find_set(ku) ) { - if( kv != idx_coc_v ) { zero_cocycles_.erase( map_it_v ); } + // Maintain the index of the 0-cocycle alive. + if (ku != idx_coc_u) { + zero_cocycles_.erase(map_it_u); + } + if (ku == dsets_.find_set(ku)) { + if (kv != idx_coc_v) { + zero_cocycles_.erase(map_it_v); + } zero_cocycles_[ku] = idx_coc_v; } } - cpx_->assign_key(sigma,cpx_->null_key()); + cpx_->assign_key(sigma, cpx_->null_key()); + } else { // If ku == kv, same connected component: create a 1-cocycle class. + create_cocycle(sigma, coeff_field_.multiplicative_identity(), coeff_field_.characteristic()); + } } - else { // If ku == kv, same connected component: create a 1-cocycle class. - create_cocycle( sigma, coeff_field_.multiplicative_identity(), coeff_field_.characteristic() ); - } -} -/* - * Compute the annotation of the boundary of a simplex. - */ -void annotation_of_the_boundary(std::map< Simplex_key, Arith_element > & map_a_ds - , Simplex_handle sigma - , int dim_sigma ) -{ - // traverses the boundary of sigma, keeps track of the annotation vectors, - // with multiplicity, in a map. - std::map < Column *, int > annotations_in_boundary; - std::pair < typename std::map< Column *, int >::iterator - , bool > result_insert_bound; - int sign = 1 - 2 * (dim_sigma % 2); // \in {-1,1} provides the sign in the - // alternate sum in the boundary. - Simplex_key key; Column * curr_col; - - for( auto sh : cpx_->boundary_simplex_range(sigma) ) - { - key = cpx_->key(sh); - if( key != cpx_->null_key() ) // A simplex with null_key is a killer, and have null annotation - { // vector. + /* + * Compute the annotation of the boundary of a simplex. + */ + void annotation_of_the_boundary( + std::map<Simplex_key, Arith_element> & map_a_ds, Simplex_handle sigma, + int dim_sigma) { + // traverses the boundary of sigma, keeps track of the annotation vectors, + // with multiplicity, in a map. + std::map<Column *, int> annotations_in_boundary; + std::pair<typename std::map<Column *, int>::iterator, bool> result_insert_bound; + int sign = 1 - 2 * (dim_sigma % 2); // \in {-1,1} provides the sign in the + // alternate sum in the boundary. + Simplex_key key; + Column * curr_col; + + for (auto sh : cpx_->boundary_simplex_range(sigma)) { + key = cpx_->key(sh); + if (key != cpx_->null_key()) { // A simplex with null_key is a killer, and have null annotation // Find its annotation vector - curr_col = ds_repr_[ dsets_.find_set(key) ]; - if( curr_col != NULL ) - { // and insert it in annotations_in_boundary with multyiplicative factor "sign". - result_insert_bound = - annotations_in_boundary.insert(std::pair<Column *,int>(curr_col,sign)); - if( !(result_insert_bound.second) ) { result_insert_bound.first->second += sign; } + curr_col = ds_repr_[dsets_.find_set(key)]; + if (curr_col != NULL) { // and insert it in annotations_in_boundary with multyiplicative factor "sign". + result_insert_bound = annotations_in_boundary.insert(std::pair<Column *, int>(curr_col, sign)); + if (!(result_insert_bound.second)) { + result_insert_bound.first->second += sign; + } } } - sign = -sign; - } - // Sum the annotations with multiplicity, using a map<key,coeff> - // to represent a sparse vector. - std::pair < typename std::map < Simplex_key, Arith_element >::iterator - , bool > result_insert_a_ds; - - for( auto ann_ref : annotations_in_boundary ) - { - if(ann_ref.second != coeff_field_.additive_identity()) // For all columns in the boundary, - { - for( auto cell_ref : ann_ref.first->col_ ) // insert every cell in map_a_ds with multiplicity - { - Arith_element w_y = - coeff_field_.times(cell_ref.coefficient_ , ann_ref.second); //coefficient * multiplicity - - if( w_y != coeff_field_.additive_identity() ) // if != 0 - { - result_insert_a_ds = map_a_ds.insert(std::pair< Simplex_key - , Arith_element >(cell_ref.key_ , w_y)); - if( !(result_insert_a_ds.second) ) //if cell_ref.key_ already a Key in map_a_ds - { - coeff_field_.plus_equal(result_insert_a_ds.first->second, w_y); - if(result_insert_a_ds.first->second == coeff_field_.additive_identity()) - { map_a_ds.erase(result_insert_a_ds.first); } + sign = -sign; + } + // Sum the annotations with multiplicity, using a map<key,coeff> + // to represent a sparse vector. + std::pair<typename std::map<Simplex_key, Arith_element>::iterator, bool> result_insert_a_ds; + + for (auto ann_ref : annotations_in_boundary) { + if (ann_ref.second != coeff_field_.additive_identity()) { // For all columns in the boundary, + for (auto cell_ref : ann_ref.first->col_) { // insert every cell in map_a_ds with multiplicity + Arith_element w_y = coeff_field_.times(cell_ref.coefficient_, ann_ref.second); // coefficient * multiplicity + + if (w_y != coeff_field_.additive_identity()) { // if != 0 + result_insert_a_ds = map_a_ds.insert(std::pair<Simplex_key, Arith_element>(cell_ref.key_, w_y)); + if (!(result_insert_a_ds.second)) { // if cell_ref.key_ already a Key in map_a_ds + result_insert_a_ds.first->second = coeff_field_.plus_equal(result_insert_a_ds.first->second, w_y); + if (result_insert_a_ds.first->second == coeff_field_.additive_identity()) { + map_a_ds.erase(result_insert_a_ds.first); + } + } } } } - } + } } -} -/* - * Update the cohomology groups under the insertion of a simplex. - */ -void update_cohomology_groups ( Simplex_handle sigma - , int dim_sigma ) -{ -//Compute the annotation of the boundary of sigma: - std::map< Simplex_key, Arith_element > map_a_ds; - annotation_of_the_boundary(map_a_ds, sigma, dim_sigma ); + /* + * Update the cohomology groups under the insertion of a simplex. + */ + void update_cohomology_groups(Simplex_handle sigma, int dim_sigma) { +// Compute the annotation of the boundary of sigma: + std::map<Simplex_key, Arith_element> map_a_ds; + annotation_of_the_boundary(map_a_ds, sigma, dim_sigma); // Update the cohomology groups: - if( map_a_ds.empty() ) { // sigma is a creator in all fields represented in coeff_field_ - if(dim_sigma < dim_max_) { create_cocycle ( sigma - , coeff_field_.multiplicative_identity() - , coeff_field_.characteristic() );} - } - else { // sigma is a destructor in at least a field in coeff_field_ - // Convert map_a_ds to a vector - A_ds_type a_ds; //admits reverse iterators - for ( auto map_a_ds_ref : map_a_ds ) - { - a_ds.push_back( std::pair< Simplex_key - , Arith_element> ( map_a_ds_ref.first - , map_a_ds_ref.second )); - } + if (map_a_ds.empty()) { // sigma is a creator in all fields represented in coeff_field_ + if (dim_sigma < dim_max_) { + create_cocycle(sigma, coeff_field_.multiplicative_identity(), + coeff_field_.characteristic()); + } + } else { // sigma is a destructor in at least a field in coeff_field_ + // Convert map_a_ds to a vector + A_ds_type a_ds; // admits reverse iterators + for (auto map_a_ds_ref : map_a_ds) { + a_ds.push_back( + std::pair<Simplex_key, Arith_element>(map_a_ds_ref.first, + map_a_ds_ref.second)); + } + + Arith_element inv_x, charac; + Arith_element prod = coeff_field_.characteristic(); // Product of characteristic of the fields + for (auto a_ds_rit = a_ds.rbegin(); + (a_ds_rit != a_ds.rend()) + && (prod != coeff_field_.multiplicative_identity()); ++a_ds_rit) { + std::tie(inv_x, charac) = coeff_field_.inverse(a_ds_rit->second, prod); - - Arith_element inv_x, charac; - Arith_element prod = coeff_field_.characteristic(); // Product of characteristic of the fields - for( auto a_ds_rit = a_ds.rbegin(); - (a_ds_rit != a_ds.rend()) && (prod != coeff_field_.multiplicative_identity()); - ++a_ds_rit ) - { - std::tie(inv_x,charac) = coeff_field_.inverse ( a_ds_rit->second - , prod ); - - if( inv_x != coeff_field_.additive_identity() ) - { - destroy_cocycle ( sigma - , a_ds - , a_ds_rit->first - , inv_x - , charac ); - prod /= charac; + if (inv_x != coeff_field_.additive_identity()) { + destroy_cocycle(sigma, a_ds, a_ds_rit->first, inv_x, charac); + prod /= charac; + } + } + if (prod != coeff_field_.multiplicative_identity() + && dim_sigma < dim_max_) { + create_cocycle(sigma, coeff_field_.multiplicative_identity(prod), prod); } } - if( prod != coeff_field_.multiplicative_identity() && dim_sigma < dim_max_ ) - { create_cocycle( sigma , coeff_field_.multiplicative_identity(prod), prod ); } } -} -/* \brief Create a new cocycle class. - * - * The class is created by the insertion of the simplex sigma. - * The methods adds a cocycle, representing the new cocycle class, - * to the matrix representing the cohomology groups. - * The new cocycle has value 0 on every simplex except on sigma - * where it worths 1.*/ -void create_cocycle ( Simplex_handle sigma - , Arith_element x - , Arith_element charac ) -{ - Simplex_key key = cpx_->key(sigma); - // Create a column containing only one cell, - Column * new_col = column_pool_->construct(Column(key)); - Cell * new_cell = cell_pool_->construct(Cell (key, x, new_col)); - new_col->col_.push_back(*new_cell); - // and insert it in the matrix, in constant time thanks to the hint cam_.end(). - // Indeed *new_col has the biggest lexicographic value because key is the - // biggest key used so far. - cam_.insert (cam_.end(), *new_col); - // Update the disjoint sets data structure. - Hcell * new_hcell = new Hcell; - new_hcell->push_back(*new_cell); - transverse_idx_[key] = cocycle(charac,new_hcell); //insert the new row - ds_repr_[key] = new_col; -} - -/* \brief Destroy a cocycle class. - * - * The cocycle class is destroyed by the insertion of sigma. - * The methods proceeds to a reduction of the matrix representing - * the cohomology groups using Gauss pivoting. The reduction zeros-out - * the row containing the cell with highest key in - * a_ds, the annotation of the boundary of simplex sigma. This key - * is "death_key".*/ -void destroy_cocycle ( Simplex_handle sigma - , A_ds_type const& a_ds - , Simplex_key death_key - , Arith_element inv_x - , Arith_element charac ) -{ - // Create a finite persistent interval - if(interval_length_policy(cpx_->simplex(death_key),sigma)) { - persistent_pairs_.push_back ( Persistent_interval ( cpx_->simplex(death_key) //creator - , sigma //destructor - , charac ) //fields - ); // for which the interval exists + /* \brief Create a new cocycle class. + * + * The class is created by the insertion of the simplex sigma. + * The methods adds a cocycle, representing the new cocycle class, + * to the matrix representing the cohomology groups. + * The new cocycle has value 0 on every simplex except on sigma + * where it worths 1.*/ + void create_cocycle(Simplex_handle sigma, Arith_element x, + Arith_element charac) { + Simplex_key key = cpx_->key(sigma); + // Create a column containing only one cell, + Column * new_col = column_pool_->construct(Column(key)); + Cell * new_cell = cell_pool_->construct(Cell(key, x, new_col)); + new_col->col_.push_back(*new_cell); + // and insert it in the matrix, in constant time thanks to the hint cam_.end(). + // Indeed *new_col has the biggest lexicographic value because key is the + // biggest key used so far. + cam_.insert(cam_.end(), *new_col); + // Update the disjoint sets data structure. + Hcell * new_hcell = new Hcell; + new_hcell->push_back(*new_cell); + transverse_idx_[key] = cocycle(charac, new_hcell); // insert the new row + ds_repr_[key] = new_col; } - auto death_key_row = transverse_idx_.find(death_key); // Find the beginning of the row. - std::pair< typename Cam::iterator, bool > result_insert_cam; - - auto row_cell_it = death_key_row->second.row_->begin(); - - while( row_cell_it != death_key_row->second.row_->end() ) // Traverse all cells in - { // the row at index death_key. - Arith_element w = coeff_field_.times_minus( inv_x , row_cell_it->coefficient_ ); - - if( w != coeff_field_.additive_identity() ) - { - Column * curr_col = row_cell_it->self_col_; ++row_cell_it; - // Disconnect the column from the rows in the CAM. - for( auto col_cell_it = curr_col->col_.begin(); - col_cell_it != curr_col->col_.end(); ++col_cell_it ) - { col_cell_it->base_hook_cam_h::unlink(); } - - // Remove the column from the CAM before modifying its value - cam_.erase( cam_.iterator_to(*curr_col) ); - // Proceed to the reduction of the column - plus_equal_column(*curr_col, a_ds, w); - - if( curr_col->col_.empty() ) // If the column is null - { - ds_repr_[ curr_col->class_key_ ] = NULL; - column_pool_->free(curr_col); //delete curr_col; - } - else - { - // Find whether the column obtained is already in the CAM - result_insert_cam = cam_.insert( *curr_col ); - if ( result_insert_cam.second ) // If it was not in the CAM before: insertion has succeeded - { - for ( auto col_cell_it = curr_col->col_.begin(); - col_cell_it != curr_col->col_.end(); ++col_cell_it ) //re-establish the row links - { transverse_idx_[ col_cell_it->key_ ].row_->push_front(*col_cell_it); } + /* \brief Destroy a cocycle class. + * + * The cocycle class is destroyed by the insertion of sigma. + * The methods proceeds to a reduction of the matrix representing + * the cohomology groups using Gauss pivoting. The reduction zeros-out + * the row containing the cell with highest key in + * a_ds, the annotation of the boundary of simplex sigma. This key + * is "death_key".*/ + void destroy_cocycle(Simplex_handle sigma, A_ds_type const& a_ds, + Simplex_key death_key, Arith_element inv_x, + Arith_element charac) { + // Create a finite persistent interval for which the interval exists + if (interval_length_policy(cpx_->simplex(death_key), sigma)) { + persistent_pairs_.push_back(Persistent_interval(cpx_->simplex(death_key) // creator + , sigma // destructor + , charac)); // fields + } + + auto death_key_row = transverse_idx_.find(death_key); // Find the beginning of the row. + std::pair<typename Cam::iterator, bool> result_insert_cam; + + auto row_cell_it = death_key_row->second.row_->begin(); + + while (row_cell_it != death_key_row->second.row_->end()) { // Traverse all cells in + // the row at index death_key. + Arith_element w = coeff_field_.times_minus(inv_x, row_cell_it->coefficient_); + + if (w != coeff_field_.additive_identity()) { + Column * curr_col = row_cell_it->self_col_; + ++row_cell_it; + // Disconnect the column from the rows in the CAM. + for (auto col_cell_it = curr_col->col_.begin(); + col_cell_it != curr_col->col_.end(); ++col_cell_it) { + col_cell_it->base_hook_cam_h::unlink(); } - else // There is already an identical column in the CAM: - { // merge two disjoint sets. - dsets_.link ( curr_col->class_key_ , - result_insert_cam.first->class_key_ ); - - Simplex_key key_tmp = dsets_.find_set( curr_col->class_key_ ); - ds_repr_[ key_tmp ] = &(*(result_insert_cam.first)); - result_insert_cam.first->class_key_ = key_tmp; - column_pool_->free(curr_col); //delete curr_col; + + // Remove the column from the CAM before modifying its value + cam_.erase(cam_.iterator_to(*curr_col)); + // Proceed to the reduction of the column + plus_equal_column(*curr_col, a_ds, w); + + if (curr_col->col_.empty()) { // If the column is null + ds_repr_[curr_col->class_key_] = NULL; + column_pool_->free(curr_col); // delete curr_col; + } else { + // Find whether the column obtained is already in the CAM + result_insert_cam = cam_.insert(*curr_col); + if (result_insert_cam.second) { // If it was not in the CAM before: insertion has succeeded + for (auto col_cell_it = curr_col->col_.begin(); col_cell_it != curr_col->col_.end(); ++col_cell_it) { + // re-establish the row links + transverse_idx_[col_cell_it->key_].row_->push_front(*col_cell_it); + } + } else { // There is already an identical column in the CAM: + // merge two disjoint sets. + dsets_.link(curr_col->class_key_, + result_insert_cam.first->class_key_); + + Simplex_key key_tmp = dsets_.find_set(curr_col->class_key_); + ds_repr_[key_tmp] = &(*(result_insert_cam.first)); + result_insert_cam.first->class_key_ = key_tmp; + column_pool_->free(curr_col); // delete curr_col; + } } - } + } else { + ++row_cell_it; + } // If w == 0, pass. } - else { ++row_cell_it; } // If w == 0, pass. - } - // Because it is a killer simplex, set the data of sigma to null_key(). - if(charac == coeff_field_.characteristic()) { cpx_->assign_key( sigma, cpx_->null_key() ); } - if(death_key_row->second.characteristics_ == charac) - { - delete death_key_row->second.row_; - transverse_idx_.erase(death_key_row); + // Because it is a killer simplex, set the data of sigma to null_key(). + if (charac == coeff_field_.characteristic()) { + cpx_->assign_key(sigma, cpx_->null_key()); + } + if (death_key_row->second.characteristics_ == charac) { + delete death_key_row->second.row_; + transverse_idx_.erase(death_key_row); + } else { + death_key_row->second.characteristics_ /= charac; + } } - else { death_key_row->second.characteristics_ /= charac; } -} -/* - * Assign: target <- target + w * other. - */ -void plus_equal_column ( Column & target - , A_ds_type const& other //value_type is pair<Simplex_key,Arith_element> - , Arith_element w ) -{ - auto target_it = target.col_.begin(); auto other_it = other.begin(); - while ( target_it != target.col_.end() && other_it != other.end() ) - { - if(target_it->key_ < other_it->first) { ++target_it; } - else { - if(target_it->key_ > other_it->first) - { - Cell * cell_tmp = cell_pool_->construct(Cell( other_it->first //key - , coeff_field_.additive_identity() - , &target)); - - coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w); - - target.col_.insert( target_it, *cell_tmp ); - - ++other_it; - } - else { //it1->key == it2->key - //target_it->coefficient_ <- target_it->coefficient_ + other_it->second * w - coeff_field_.plus_times_equal( target_it->coefficient_ , other_it->second , w); - if( target_it->coefficient_ == coeff_field_.additive_identity() ) - { - auto tmp_it = target_it; - ++target_it; ++other_it; // iterators remain valid - Cell * tmp_cell_ptr = &(*tmp_it); - target.col_.erase(tmp_it); // removed from column - - coeff_field_.clear_coefficient(tmp_cell_ptr->coefficient_); - cell_pool_->free(tmp_cell_ptr); // delete from memory + /* + * Assign: target <- target + w * other. + */ + void plus_equal_column(Column & target, A_ds_type const& other // value_type is pair<Simplex_key,Arith_element> + , Arith_element w) { + auto target_it = target.col_.begin(); + auto other_it = other.begin(); + while (target_it != target.col_.end() && other_it != other.end()) { + if (target_it->key_ < other_it->first) { + ++target_it; + } else { + if (target_it->key_ > other_it->first) { + Cell * cell_tmp = cell_pool_->construct(Cell(other_it->first // key + , coeff_field_.additive_identity(), &target)); + + cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w); + + target.col_.insert(target_it, *cell_tmp); + + ++other_it; + } else { // it1->key == it2->key + // target_it->coefficient_ <- target_it->coefficient_ + other_it->second * w + target_it->coefficient_ = coeff_field_.plus_times_equal(target_it->coefficient_, other_it->second, w); + if (target_it->coefficient_ == coeff_field_.additive_identity()) { + auto tmp_it = target_it; + ++target_it; + ++other_it; // iterators remain valid + Cell * tmp_cell_ptr = &(*tmp_it); + target.col_.erase(tmp_it); // removed from column + + coeff_field_.clear_coefficient(tmp_cell_ptr->coefficient_); + cell_pool_->free(tmp_cell_ptr); // delete from memory + } else { + ++target_it; + ++other_it; + } } - else { ++target_it; ++other_it; } } } - } - while(other_it != other.end()) - { - Cell * cell_tmp = cell_pool_->construct(Cell( other_it->first //key - , coeff_field_.additive_identity() - , &target)); + while (other_it != other.end()) { + Cell * cell_tmp = cell_pool_->construct(Cell(other_it->first, coeff_field_.additive_identity(), &target)); + cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w); + target.col_.insert(target.col_.end(), *cell_tmp); - coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w); + ++other_it; + } + } - target.col_.insert( target.col_.end(), *cell_tmp ); + /* + * Compare two intervals by length. + */ + struct cmp_intervals_by_length { + explicit cmp_intervals_by_length(Complex_ds * sc) + : sc_(sc) { + } + bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) { + return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1)) + > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2))); + } + Complex_ds * sc_; + }; - ++other_it; + public: + /** \brief Output the persistence diagram in ostream. + * + * The file format is the following: + * p1*...*pr dim b d + * + * where "dim" is the dimension of the homological feature, + * b and d are respectively the birth and death of the feature and + * p1*...*pr is the product of prime numbers pi such that the homology + * feature exists in homology with Z/piZ coefficients. + */ + void output_diagram(std::ostream& ostream = std::cout) { + cmp_intervals_by_length cmp(cpx_); + persistent_pairs_.sort(cmp); + for (auto pair : persistent_pairs_) { + ostream << get<2>(pair) << " " << cpx_->dimension(get<0>(pair)) << " " + << cpx_->filtration(get<0>(pair)) << " " + << cpx_->filtration(get<1>(pair)) << " " << std::endl; + } } -} -/* - * Compare two intervals by length. - */ -struct cmp_intervals_by_length { - cmp_intervals_by_length( Complex_ds * sc ) : sc_ (sc) {} - bool operator() ( Persistent_interval & p1 - , Persistent_interval & p2 ) - { - return ( sc_->filtration( get<1>(p1) ) - sc_->filtration( get<0>(p1) ) - > sc_->filtration( get<1>(p2) ) - sc_->filtration( get<0>(p2) ) ); - } - Complex_ds * sc_; -}; + private: + /* + * Structure representing a cocycle. + */ + struct cocycle { + cocycle() { + } + cocycle(Arith_element characteristics, Hcell * row) + : row_(row), + characteristics_(characteristics) { + } -public: -/** \brief Output the persistence diagram in ostream. - * - * The file format is the following: - * p1*...*pr dim b d - * - * where "dim" is the dimension of the homological feature, - * b and d are respectively the birth and death of the feature and - * p1*...*pr is the product of prime numbers pi such that the homology - * feature exists in homology with Z/piZ coefficients. - */ -void output_diagram(std::ostream& ostream = std::cout) -{ - cmp_intervals_by_length cmp( cpx_ ); - persistent_pairs_.sort( cmp ); - for(auto pair : persistent_pairs_) - { - ostream << get<2>(pair) << " " - << cpx_->dimension(get<0>(pair)) << " " - << cpx_->filtration(get<0>(pair)) << " " - << cpx_->filtration(get<1>(pair)) << " " - << std::endl; - } -} + Hcell * row_; // points to the corresponding row in the CAM + Arith_element characteristics_; // product of field characteristics for which the cocycle exist + }; -private: -/* - * Structure representing a cocycle. - */ -struct cocycle { - cocycle() {} - cocycle( Arith_element characteristics - , Hcell * row ) - : row_(row), characteristics_(characteristics) {} - - Hcell * row_; //points to the corresponding row in the CAM - Arith_element characteristics_; //product of field characteristics for which the cocycle exist + public: + Complex_ds * cpx_; + int dim_max_; + CoefficientField coeff_field_; + + /* Disjoint sets data structure to link the model of FilteredComplex + * with the compressed annotation matrix. + * ds_rank_ is a property map Simplex_key -> int, ds_parent_ is a property map + * Simplex_key -> simplex_key_t */ + std::vector<int> ds_rank_; + std::vector<Simplex_key> ds_parent_; + std::vector<Column *> ds_repr_; + boost::disjoint_sets<int *, Simplex_key *> dsets_; + /* The compressed annotation matrix fields.*/ + Cam cam_; + /* Dictionary establishing the correspondance between the Simplex_key of + * the root vertex in the union-find ds and the Simplex_key of the vertex which + * created the connected component as a 0-dimension homology feature.*/ + std::map<Simplex_key, Simplex_key> zero_cocycles_; + /* Key -> row. */ + std::map<Simplex_key, cocycle> transverse_idx_; + /* Persistent intervals. */ + std::list<Persistent_interval> persistent_pairs_; + length_interval interval_length_policy; + + boost::object_pool<Column> * column_pool_; + boost::object_pool<Cell> * cell_pool_; }; -public: - Complex_ds * cpx_; - int dim_max_; - CoefficientField coeff_field_; - -/* Disjoint sets data structure to link the model of FilteredComplex - * with the compressed annotation matrix. - * ds_rank_ is a property map Simplex_key -> int, ds_parent_ is a property map - * Simplex_key -> simplex_key_t */ - std::vector< int > ds_rank_; - std::vector< Simplex_key > ds_parent_; - std::vector< Column * > ds_repr_; - boost::disjoint_sets< int *, Simplex_key * > dsets_; -/* The compressed annotation matrix fields.*/ - Cam cam_; -/* Dictionary establishing the correspondance between the Simplex_key of - * the root vertex in the union-find ds and the Simplex_key of the vertex which - * created the connected component as a 0-dimension homology feature.*/ - std::map<Simplex_key,Simplex_key> zero_cocycles_; -/* Key -> row. */ - std::map< Simplex_key , cocycle > transverse_idx_; -/* Persistent intervals. */ - std::list< Persistent_interval > persistent_pairs_; - length_interval interval_length_policy; - - boost::object_pool< Column > * column_pool_; - boost::object_pool< Cell > * cell_pool_; -}; +/** @} */ // end defgroup persistent_cohomology -/** @} */ //end defgroup persistent_cohomology +} // namespace persistent_cohomology -} // namespace GUDHI +} // namespace Gudhi -#endif // _PERSISTENCECOMPUTATION_SIMPLEXTREE_ +#endif // SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_H_ diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h index af0d6605..419bd2eb 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h @@ -1,106 +1,123 @@ - /* 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): Clément Maria - * - * 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 GUDHI_FIELD_ZP_H -#define GUDHI_FIELD_ZP_H - -namespace Gudhi{ +/* 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): Clément Maria + * + * 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 SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_FIELD_ZP_H_ +#define SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_FIELD_ZP_H_ + +#include <utility> +#include <vector> + +namespace Gudhi { + +namespace persistent_cohomology { /** \brief Structure representing the coefficient field \f$\mathbb{Z}/p\mathbb{Z}\f$ - * - * \implements CoefficientField - * \ingroup persistent_cohomology - */ + * + * \implements CoefficientField + * \ingroup persistent_cohomology + */ class Field_Zp { -public: -typedef int Element; - -Field_Zp() -: Prime(-1) -, inverse_() {} - -void init(int charac ) { - assert(charac <= 32768); - Prime = charac; - inverse_.clear(); - inverse_.reserve(charac); - inverse_.push_back(0); - for(int i=1 ; i<Prime ; ++i) - { - int inv = 1; - while(((inv * i) % Prime) != 1) ++inv; - inverse_.push_back(inv); + public: + typedef int Element; + + Field_Zp() + : Prime(0), + inverse_(), + mult_id_all(1), + add_id_all(0) { } -} -/** Set x <- x + w * y*/ -void plus_times_equal ( Element & x, Element y, Element w ) -{ x = (x + w * y) % Prime; } - -// operator= defined on Element - -/** Returns y * w */ -Element times ( Element y, int w ) { - Element res = (y * w) % Prime; - if(res < 0) return res+Prime; - else return res; -} - -void clear_coefficient(Element x) {} + void init(uint16_t charac) { + assert(charac != 0); // division by zero + Prime = charac; + inverse_.clear(); + inverse_.reserve(charac); + inverse_.push_back(0); + for (int i = 1; i < Prime; ++i) { + int inv = 1; + while (((inv * i) % Prime) != 1) + ++inv; + inverse_.push_back(inv); + } + } -void plus_equal(Element & x, Element y) { x = ((x+y)%Prime); } + /** Set x <- x + w * y*/ + Element plus_times_equal(const Element& x, const Element& y, const Element& w) { + assert(Prime != 0); // division by zero + Element result = (x + w * y) % Prime; + if (result < 0) + result += Prime; + return result; + } -/** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ -Element additive_identity () { return 0; } -/** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ -Element multiplicative_identity ( Element P = 0) { return 1; } -/** Returns the inverse in the field. Modifies P.*/ -std::pair<Element,Element> inverse ( Element x - , Element P ) -{ return std::pair<Element,Element>(inverse_[x],P); -} // <------ return the product of field characteristic for which x is invertible +// operator= defined on Element -/** Returns -x * y.*/ -Element times_minus ( Element x, Element y ) -{ - Element out = (-x * y) % Prime; - return (out < 0) ? out + Prime : out; -} + /** Returns y * w */ + Element times(const Element& y, const Element& w) { + return plus_times_equal(0, y, (Element)w); + } + void clear_coefficient(Element x) { + } -bool is_one ( Element x ) { return x == 1; } -bool is_zero ( Element x ) { return x == 0; } + Element plus_equal(const Element& x, const Element& y) { + return plus_times_equal(x, y, (Element)1); + } -//bool is_null() + /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ + const Element& additive_identity() const { + return add_id_all; + } + /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ + const Element& multiplicative_identity(Element P = 0) const { + return mult_id_all; + } + /** Returns the inverse in the field. Modifies P.*/ + std::pair<Element, Element> inverse(Element x, Element P) { + return std::pair<Element, Element>(inverse_[x], P); + } // <------ return the product of field characteristic for which x is invertible + + /** Returns -x * y.*/ + Element times_minus(Element x, Element y) { + assert(Prime != 0); // division by zero + Element out = (-x * y) % Prime; + return (out < 0) ? out + Prime : out; + } -/** \brief Returns the characteristic \f$p\f$ of the field.*/ -Element characteristic() { return Prime; } + /** \brief Returns the characteristic \f$p\f$ of the field.*/ + const uint16_t& characteristic() const { + return Prime; + } -private: - Element Prime; -/** Property map Element -> Element, which associate to an element its inverse in the field.*/ - std::vector< Element > inverse_; + private: + uint16_t Prime; + /** Property map Element -> Element, which associate to an element its inverse in the field.*/ + std::vector<Element> inverse_; + const Element mult_id_all; + const Element add_id_all; }; -} // namespace GUDHI +} // namespace persistent_cohomology + +} // namespace Gudhi -#endif // GUDHI_FIELD_ZP_H +#endif // SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_FIELD_ZP_H_ diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h index 9dd0998c..91937c65 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h @@ -1,167 +1,189 @@ - /* 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): Clément Maria - * - * 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 GUDHI_MULTI_FIELD_H -#define GUDHI_MULTI_FIELD_H - -#include <iostream> -#include <vector> +/* 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): Clément Maria + * + * 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 SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_ +#define SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_ + #include <gmpxx.h> -namespace Gudhi{ +#include <vector> +#include <utility> + +namespace Gudhi { + +namespace persistent_cohomology { /** \brief Structure representing coefficients in a set of finite fields simultaneously - * using the chinese remainder theorem. - * - * \implements CoefficientField - * \ingroup persistent_cohomology - - * Details on the algorithms may be found in \cite boissonnat:hal-00922572 - */ + * using the chinese remainder theorem. + * + * \implements CoefficientField + * \ingroup persistent_cohomology + + * Details on the algorithms may be found in \cite boissonnat:hal-00922572 + */ class Multi_field { -public: - typedef mpz_class Element; - - Multi_field () {} - -/* Initialize the multi-field. The generation of prime numbers might fail with - * a very small probability.*/ - void init(int min_prime, int max_prime) - { - if(max_prime<2) - { std::cerr << "There is no prime less than " << max_prime << std::endl; } - if(min_prime > max_prime) - { std::cerr << "No prime in ["<<min_prime<<":"<<max_prime<<"]"<<std::endl; } + public: + typedef mpz_class Element; + + Multi_field() + : prod_characteristics_(0), + mult_id_all(0), + add_id_all(0) { + } + + /* Initialize the multi-field. The generation of prime numbers might fail with + * a very small probability.*/ + void init(uint16_t min_prime, uint16_t max_prime) { + if (max_prime < 2) { + std::cerr << "There is no prime less than " << max_prime << std::endl; + } + if (min_prime > max_prime) { + std::cerr << "No prime in [" << min_prime << ":" << max_prime << "]" + << std::endl; + } // fill the list of prime numbers - unsigned int curr_prime = min_prime; - mpz_t tmp_prime; mpz_init_set_ui(tmp_prime,min_prime); - //test if min_prime is prime - int is_prime = mpz_probab_prime_p(tmp_prime,25); //probabilistic primality test - - if(is_prime == 0) //min_prime is composite - { - mpz_nextprime(tmp_prime,tmp_prime); + uint16_t curr_prime = min_prime; + mpz_t tmp_prime; + mpz_init_set_ui(tmp_prime, min_prime); + // test if min_prime is prime + int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test + + if (is_prime == 0) { // min_prime is composite + mpz_nextprime(tmp_prime, tmp_prime); curr_prime = mpz_get_ui(tmp_prime); } - - while (curr_prime <= max_prime) - { + + while (curr_prime <= max_prime) { primes_.push_back(curr_prime); - mpz_nextprime(tmp_prime,tmp_prime); + mpz_nextprime(tmp_prime, tmp_prime); curr_prime = mpz_get_ui(tmp_prime); } - //set m to primorial(bound_prime) + // set m to primorial(bound_prime) prod_characteristics_ = 1; - for(auto p : primes_) - { mpz_mul_ui(prod_characteristics_.get_mpz_t(), - prod_characteristics_.get_mpz_t(), - p); + for (auto p : primes_) { + mpz_mul_ui(prod_characteristics_.get_mpz_t(), + prod_characteristics_.get_mpz_t(), p); } - num_primes_ = primes_.size(); - - //Uvect_ - Element Ui; Element tmp_elem; - for(auto p : primes_) - { + // Uvect_ + Element Ui; + Element tmp_elem; + for (auto p : primes_) { + assert(p != 0); // division by zero tmp_elem = prod_characteristics_ / p; - //Element tmp_elem_bis = 10; - mpz_powm_ui ( tmp_elem.get_mpz_t() - , tmp_elem.get_mpz_t() - , p - 1 - , prod_characteristics_.get_mpz_t() ); + // Element tmp_elem_bis = 10; + mpz_powm_ui(tmp_elem.get_mpz_t(), tmp_elem.get_mpz_t(), p - 1, + prod_characteristics_.get_mpz_t()); Uvect_.push_back(tmp_elem); } mult_id_all = 0; - for(int idx = 0; idx < num_primes_; ++idx) - { mult_id_all = (mult_id_all + Uvect_[idx]) % prod_characteristics_; } - + for (auto uvect : Uvect_) { + assert(prod_characteristics_ != 0); // division by zero + mult_id_all = (mult_id_all + uvect) % prod_characteristics_; + } } - void clear_coefficient(Element & x) { mpz_clear(x.get_mpz_t()); } + void clear_coefficient(Element & x) { + mpz_clear(x.get_mpz_t()); + } /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ - Element additive_identity () { return 0; } + const Element& additive_identity() const { + return add_id_all; + } /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ - Element multiplicative_identity () { return mult_id_all; }// 1 everywhere + const Element& multiplicative_identity() const { + return mult_id_all; + } // 1 everywhere - Element multiplicative_identity (Element Q) - { - if(Q == prod_characteristics_) { return multiplicative_identity(); } + Element multiplicative_identity(Element Q) { + if (Q == prod_characteristics_) { + return multiplicative_identity(); + } + assert(prod_characteristics_ != 0); // division by zero Element mult_id = 0; - for(int idx = 0; idx < num_primes_; ++idx) { - if( (Q % primes_[idx]) == 0 ) - { mult_id = (mult_id + Uvect_[idx]) % prod_characteristics_; } + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + assert(primes_[idx] != 0); // division by zero + if ((Q % primes_[idx]) == 0) { + mult_id = (mult_id + Uvect_[idx]) % prod_characteristics_; + } } return mult_id; - } + } /** Returns y * w */ - Element times ( Element y, int w ) { - Element tmp = (y*w) % prod_characteristics_; - if(tmp < 0) return prod_characteristics_ + tmp; - return tmp; + Element times(const Element& y, const Element& w) { + return plus_times_equal(0, y, w); } - void plus_equal(Element & x, Element y) - { x += y; x %= prod_characteristics_; } + Element plus_equal(const Element& x, const Element& y) { + return plus_times_equal(x, y, (Element)1); + } /** \brief Returns the characteristic \f$p\f$ of the field.*/ - Element characteristic() { return prod_characteristics_; } + const Element& characteristic() const { + return prod_characteristics_; + } /** Returns the inverse in the field. Modifies P.*/ - std::pair<Element,Element> inverse ( Element x - , Element QS ) - { + std::pair<Element, Element> inverse(Element x, Element QS) { Element QR; - mpz_gcd( QR.get_mpz_t(), x.get_mpz_t(), QS.get_mpz_t() ); // QR <- gcd(x,QS) - if( QR == QS ) return std::pair<Element,Element>(additive_identity() - , multiplicative_identity() ); //partial inverse is 0 + mpz_gcd(QR.get_mpz_t(), x.get_mpz_t(), QS.get_mpz_t()); // QR <- gcd(x,QS) + if (QR == QS) + return std::pair<Element, Element>(additive_identity(), multiplicative_identity()); // partial inverse is 0 Element QT = QS / QR; Element inv_qt; mpz_invert(inv_qt.get_mpz_t(), x.get_mpz_t(), QT.get_mpz_t()); - return std::pair<Element,Element>( - (inv_qt * multiplicative_identity(QT)) % prod_characteristics_ - , QT ); + assert(prod_characteristics_ != 0); // division by zero + return std::pair<Element, Element>( + (inv_qt * multiplicative_identity(QT)) % prod_characteristics_, QT); } /** Returns -x * y.*/ - Element times_minus ( Element x, Element y ) - { return prod_characteristics_ - ((x*y)%prod_characteristics_); } + Element times_minus(const Element& x, const Element& y) { + assert(prod_characteristics_ != 0); // division by zero + return prod_characteristics_ - ((x * y) % prod_characteristics_); + } /** Set x <- x + w * y*/ - void plus_times_equal ( Element & x, Element y, Element w ) - { x = (x + w * y) % prod_characteristics_; } - - Element prod_characteristics_; //product of characteristics of the fields - //represented by the multi-field class - std::vector<int> primes_; //all the characteristics of the fields - std::vector<Element> Uvect_; - size_t num_primes_; //number of fields - Element mult_id_all; + Element plus_times_equal(const Element& x, const Element& y, const Element& w) { + assert(prod_characteristics_ != 0); // division by zero + Element result = (x + w * y) % prod_characteristics_; + if (result < 0) + result += prod_characteristics_; + return result; + } + Element prod_characteristics_; // product of characteristics of the fields + // represented by the multi-field class + std::vector<uint16_t> primes_; // all the characteristics of the fields + std::vector<Element> Uvect_; + Element mult_id_all; + const Element add_id_all; }; -} // namespace GUDHI +} // namespace persistent_cohomology + +} // namespace Gudhi -#endif // GUDHI_MULTI_FIELD_H +#endif // SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_ diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h index 0702c58e..fcec819a 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h @@ -1,153 +1,148 @@ - /* 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): Clément Maria - * - * 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 GUDHI_COLUMN_LIST_H -#define GUDHI_COLUMN_LIST_H +/* 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): Clément Maria + * + * 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 SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ +#define SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ + +#include <list> #include "boost/tuple/tuple.hpp" #include "boost/intrusive/set.hpp" #include "boost/intrusive/list.hpp" -namespace Gudhi{ +namespace Gudhi { -template < typename SimplexKey - , typename ArithmeticElement - > +namespace persistent_cohomology { + +template<typename SimplexKey, typename ArithmeticElement> class Persistent_cohomology_column; -struct cam_h_tag; // for horizontal traversal in the CAM -struct cam_v_tag; // for vertical traversal in the CAM +struct cam_h_tag; +// for horizontal traversal in the CAM +struct cam_v_tag; +// for vertical traversal in the CAM -typedef boost::intrusive::list_base_hook - < boost::intrusive::tag < cam_h_tag > - , boost::intrusive::link_mode < boost::intrusive::auto_unlink > //allows .unlink() - > base_hook_cam_h; +typedef boost::intrusive::list_base_hook<boost::intrusive::tag<cam_h_tag>, + boost::intrusive::link_mode<boost::intrusive::auto_unlink> // allows .unlink() +> base_hook_cam_h; -typedef boost::intrusive::list_base_hook - < boost::intrusive::tag < cam_v_tag > - , boost::intrusive::link_mode < boost::intrusive::normal_link > //faster hook, less safe - > base_hook_cam_v; +typedef boost::intrusive::list_base_hook<boost::intrusive::tag<cam_v_tag>, + boost::intrusive::link_mode<boost::intrusive::normal_link> // faster hook, less safe +> base_hook_cam_v; /** \internal - * \brief - * - */ -template < typename SimplexKey - , typename ArithmeticElement - > -class Persistent_cohomology_cell -: public base_hook_cam_h -, public base_hook_cam_v -{ - public: - template < class T1, class T2 > friend class Persistent_cohomology; - friend class Persistent_cohomology_column < SimplexKey , ArithmeticElement >; - - typedef Persistent_cohomology_column< SimplexKey, ArithmeticElement > Column; - - Persistent_cohomology_cell( SimplexKey key - , ArithmeticElement x - , Column * self_col) - : key_(key) - , coefficient_(x) - , self_col_(self_col) {} - - SimplexKey key_; - ArithmeticElement coefficient_; - Column * self_col_; -}; - - - + * \brief + * + */ +template<typename SimplexKey, typename ArithmeticElement> +class Persistent_cohomology_cell : public base_hook_cam_h, + public base_hook_cam_v { + public: + template<class T1, class T2> friend class Persistent_cohomology; + friend class Persistent_cohomology_column<SimplexKey, ArithmeticElement>; + + typedef Persistent_cohomology_column<SimplexKey, ArithmeticElement> Column; + + Persistent_cohomology_cell(SimplexKey key, ArithmeticElement x, + Column * self_col) + : key_(key), + coefficient_(x), + self_col_(self_col) { + } + SimplexKey key_; + ArithmeticElement coefficient_; + Column * self_col_; +}; /* - * \brief Sparse column for the Compressed Annotation Matrix. - * - * The non-zero coefficients of the column are stored in a - * boost::intrusive::list. Contains a hook to be stored in a - * boost::intrusive::set. - */ -template < typename SimplexKey - , typename ArithmeticElement > -class Persistent_cohomology_column -: public boost::intrusive::set_base_hook - < boost::intrusive::link_mode< boost::intrusive::normal_link > > -{ -private: - template < class T1, class T2 > friend class Persistent_cohomology; - - typedef Persistent_cohomology_cell < SimplexKey, ArithmeticElement > Cell; - typedef boost::intrusive::list < Cell - , boost::intrusive::constant_time_size<false> - , boost::intrusive::base_hook< base_hook_cam_v > - > Col_type; - -/** \brief Creates an empty column.*/ - Persistent_cohomology_column (SimplexKey key) - { + * \brief Sparse column for the Compressed Annotation Matrix. + * + * The non-zero coefficients of the column are stored in a + * boost::intrusive::list. Contains a hook to be stored in a + * boost::intrusive::set. + */ +template<typename SimplexKey, typename ArithmeticElement> +class Persistent_cohomology_column : public boost::intrusive::set_base_hook< + boost::intrusive::link_mode<boost::intrusive::normal_link> > { + private: + template<class T1, class T2> friend class Persistent_cohomology; + + typedef Persistent_cohomology_cell<SimplexKey, ArithmeticElement> Cell; + typedef boost::intrusive::list<Cell, + boost::intrusive::constant_time_size<false>, + boost::intrusive::base_hook<base_hook_cam_v> > Col_type; + + /** \brief Creates an empty column.*/ + explicit Persistent_cohomology_column(SimplexKey key) { class_key_ = key; col_ = Col_type(); } -public: - /** Copy constructor.*/ - Persistent_cohomology_column( Persistent_cohomology_column const &other ) - : col_() - , class_key_(other.class_key_) - { if(!other.col_.empty()) std::cerr << "Copying a non-empty column.\n"; } - -/** \brief Returns true iff the column is null.*/ - bool is_null() { return col_.empty(); } -/** \brief Returns the key of the representative simplex of - * the set of simplices having this column as annotation vector - * in the compressed annotation matrix.*/ - SimplexKey class_key () { return class_key_; } - -/** \brief Lexicographic comparison of two columns.*/ -friend bool operator< ( const Persistent_cohomology_column& c1 - , const Persistent_cohomology_column& c2) - { - typename Col_type::const_iterator it1 = c1.col_.begin(); - typename Col_type::const_iterator it2 = c2.col_.begin(); - while(it1 != c1.col_.end() && it2 != c2.col_.end()) - { - if(it1->key_ == it2->key_) - { if(it1->coefficient_ == it2->coefficient_) { ++it1; ++it2; } - else { return it1->coefficient_ < it2->coefficient_; } } - else { return it1->key_ < it2->key_; } - } - return (it2 != c2.col_.end()); + public: + /** Copy constructor.*/ + Persistent_cohomology_column(Persistent_cohomology_column const &other) + : col_(), + class_key_(other.class_key_) { + if (!other.col_.empty()) + std::cerr << "Copying a non-empty column.\n"; + } + + /** \brief Returns true iff the column is null.*/ + bool is_null() { + return col_.empty(); + } + /** \brief Returns the key of the representative simplex of + * the set of simplices having this column as annotation vector + * in the compressed annotation matrix.*/ + SimplexKey class_key() { + return class_key_; } - // void display() - // { - // for(auto cell : col_) - // { std::cout << "(" << cell.key_ <<":"<<cell.coefficient_<<") "; } - // } + /** \brief Lexicographic comparison of two columns.*/ + friend bool operator<(const Persistent_cohomology_column& c1, + const Persistent_cohomology_column& c2) { + typename Col_type::const_iterator it1 = c1.col_.begin(); + typename Col_type::const_iterator it2 = c2.col_.begin(); + while (it1 != c1.col_.end() && it2 != c2.col_.end()) { + if (it1->key_ == it2->key_) { + if (it1->coefficient_ == it2->coefficient_) { + ++it1; + ++it2; + } else { + return it1->coefficient_ < it2->coefficient_; + } + } else { + return it1->key_ < it2->key_; + } + } + return (it2 != c2.col_.end()); + } - Col_type col_; - SimplexKey class_key_; + Col_type col_; + SimplexKey class_key_; }; -} // namespace GUDHI +} // namespace persistent_cohomology + +} // namespace Gudhi -#endif // GUDHI_COLUMN_LIST_H +#endif // SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ diff --git a/src/Persistent_cohomology/test/CMakeLists.txt b/src/Persistent_cohomology/test/CMakeLists.txt new file mode 100644 index 00000000..80689113 --- /dev/null +++ b/src/Persistent_cohomology/test/CMakeLists.txt @@ -0,0 +1,29 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHITestSimplexTree) + +if(NOT MSVC) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +endif() + +add_executable ( persistent_cohomology_unit_test persistent_cohomology_unit_test.cpp ) +target_link_libraries(persistent_cohomology_unit_test ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +# Unitary tests +add_test(persistent_cohomology_unit_test ${CMAKE_CURRENT_BINARY_DIR}/persistent_cohomology_unit_test) + +if(GMPXX_FOUND AND GMP_FOUND) + add_executable ( persistent_cohomology_unit_test_multi_field persistent_cohomology_unit_test_multi_field.cpp ) + target_link_libraries(persistent_cohomology_unit_test_multi_field ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + + # Unitary tests + #add_test(persistent_cohomology_unit_test_multi_field ${CMAKE_CURRENT_BINARY_DIR}/persistent_cohomology_unit_test_multi_field) +endif() + +if (LCOV_PATH) + # Lcov code coverage of unitary test + add_test(src/Persistent_cohomology/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Persistent_cohomology) +endif() + +cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Persistent_cohomology/include/gudhi") diff --git a/src/Persistent_cohomology/test/README b/src/Persistent_cohomology/test/README new file mode 100644 index 00000000..56474fa0 --- /dev/null +++ b/src/Persistent_cohomology/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./persistent_cohomology_unit_test --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp new file mode 100644 index 00000000..1b0cc152 --- /dev/null +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp @@ -0,0 +1,171 @@ +#define BOOST_TEST_MODULE const_string test +#include <boost/test/included/unit_test.hpp> +#include <boost/system/error_code.hpp> +#include <boost/chrono/thread_clock.hpp> +#include <iostream> +#include <string> + +#include <utility> // std::pair, std::make_pair + +#include <cmath> // float comparison +#include <limits> + +#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/reader_utils.h" +#include "gudhi/Simplex_tree.h" +#include "gudhi/Persistent_cohomology.h" + +using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; + +typedef Simplex_tree<> typeST; + +std::string test_rips_persistence(int coefficient, int min_persistence) { + const std::string inputFile("simplex_tree_file_for_unit_test.txt"); + std::ifstream simplex_tree_stream; + simplex_tree_stream.open(inputFile.c_str()); + typeST st; + simplex_tree_stream >> st; + simplex_tree_stream.close(); + + // Display the Simplex_tree + std::cout << "The complex contains " << st.num_simplices() << " simplices" << " - dimension= " << st.dimension() + << " - filtration= " << st.filtration() << std::endl; + + // Check + BOOST_CHECK(st.num_simplices() == 98); + BOOST_CHECK(st.dimension() == 3); + BOOST_CHECK(st.filtration() == 1.89); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology<Simplex_tree<>, Field_Zp> pcoh(st); + + pcoh.init_coefficients( coefficient ); // initializes the coefficient field for homology + // Check infinite rips + pcoh.compute_persistent_cohomology( min_persistence ); // Minimal lifetime of homology feature to be recorded. + std::ostringstream ossInfinite; + + pcoh.output_diagram(ossInfinite); + std::string strInfinite = ossInfinite.str(); + return strInfinite; +} + +void test_rips_persistence_in_dimension(int dimension) { + std::string value0(" 0 0.02 1.12"); + std::string value1(" 0 0.03 1.13"); + std::string value2(" 0 0.04 1.14"); + std::string value3(" 0 0.05 1.15"); + std::string value4(" 0 0.06 1.16"); + std::string value5(" 0 0.07 1.17"); + std::string value6(" 0 0.08 1.18"); + std::string value7(" 0 0.09 1.19"); + std::string value8(" 0 0 inf" ); + std::string value9(" 0 0.01 inf" ); + + value0.insert(0,std::to_string(dimension)); + value1.insert(0,std::to_string(dimension)); + value2.insert(0,std::to_string(dimension)); + value3.insert(0,std::to_string(dimension)); + value4.insert(0,std::to_string(dimension)); + value5.insert(0,std::to_string(dimension)); + value6.insert(0,std::to_string(dimension)); + value7.insert(0,std::to_string(dimension)); + value8.insert(0,std::to_string(dimension)); + value9.insert(0,std::to_string(dimension)); + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_SINGLE_FIELD DIM=" << dimension << " MIN_PERS=0" << std::endl; + + std::string str_rips_persistence = test_rips_persistence(dimension, 0); + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value4) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value5) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value6) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value7) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value8) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value9) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_SINGLE_FIELD DIM=" << dimension << " MIN_PERS=1" << std::endl; + + str_rips_persistence = test_rips_persistence(dimension, 1); + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value4) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value5) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value6) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value7) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value8) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value9) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_SINGLE_FIELD DIM=" << dimension << " MIN_PERS=2" << std::endl; + + str_rips_persistence = test_rips_persistence(dimension, 2); + + BOOST_CHECK(str_rips_persistence.find(value0) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value1) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value2) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value3) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value4) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value5) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value6) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value7) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value8) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value9) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_SINGLE_FIELD DIM=" << dimension << " MIN_PERS=Inf" << std::endl; + + str_rips_persistence = test_rips_persistence(dimension, std::numeric_limits<int>::max()); + + BOOST_CHECK(str_rips_persistence.find(value0) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value1) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value2) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value3) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value4) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value5) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value6) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value7) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value8) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value9) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_1 ) +{ + test_rips_persistence_in_dimension(1); +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_2 ) +{ + test_rips_persistence_in_dimension(2); +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_3 ) +{ + test_rips_persistence_in_dimension(3); +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_5 ) +{ + test_rips_persistence_in_dimension(5); +} + +// TODO(VR): not working from 6 +// std::string str_rips_persistence = test_rips_persistence(6, 0); +// TODO(VR): division by zero +// std::string str_rips_persistence = test_rips_persistence(0, 0); diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp new file mode 100644 index 00000000..91cf5784 --- /dev/null +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp @@ -0,0 +1,130 @@ +#define BOOST_TEST_MODULE const_string test +#include <boost/test/included/unit_test.hpp> +#include <boost/system/error_code.hpp> +#include <boost/chrono/thread_clock.hpp> +#include <iostream> +#include <string> + +#include <utility> // std::pair, std::make_pair + +#include <cmath> // float comparison +#include <limits> + +#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/reader_utils.h" +#include "gudhi/Simplex_tree.h" +#include "gudhi/Persistent_cohomology.h" +#include "gudhi/Persistent_cohomology/Multi_field.h" + +using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; + +typedef Simplex_tree<> typeST; + +std::string test_rips_persistence(int min_coefficient, int max_coefficient, int min_persistence) { + const std::string inputFile("simplex_tree_file_for_multi_field_unit_test.txt"); + std::ifstream simplex_tree_stream; + simplex_tree_stream.open(inputFile.c_str()); + typeST st; + simplex_tree_stream >> st; + simplex_tree_stream.close(); + + // Display the Simplex_tree + std::cout << "The complex contains " << st.num_simplices() << " simplices" << " - dimension= " << st.dimension() + << " - filtration= " << st.filtration() << std::endl; + + // Check + BOOST_CHECK(st.num_simplices() == 6142604); + BOOST_CHECK(st.dimension() == 3); + BOOST_CHECK(st.filtration() == 0.249999); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology<Simplex_tree<>, Multi_field> pcoh(st); + + pcoh.init_coefficients(min_coefficient, max_coefficient); // initializes the coefficient field for homology + // Check infinite rips + pcoh.compute_persistent_cohomology(min_persistence); // Minimal lifetime of homology feature to be recorded. + + std::ostringstream ossRips; + pcoh.output_diagram(ossRips); + + std::string strRips = ossRips.str(); + return strRips; +} + +void test_rips_persistence_in_dimension(int min_dimension, int max_dimension) { + std::string value0(" 0 0 inf"); + std::string value1(" 1 0.0702103 inf"); + std::string value2("2 1 0.0702103 inf"); + std::string value3("2 2 0.159992 inf"); + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD MIN_DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=0" << std::endl; + + std::string str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, 1); + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=2" << std::endl; + + str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, 2); + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=3" << std::endl; + + str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, 3); + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=Inf" << std::endl; + + str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, std::numeric_limits<int>::max()); + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_multi_field_dim_1_2 ) +{ + test_rips_persistence_in_dimension(1, 2); +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_multi_field_dim_2_3 ) +{ + test_rips_persistence_in_dimension(2, 3); +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_multi_field_dim_1_5 ) +{ + test_rips_persistence_in_dimension(1, 5); +} + +// TODO(VR): not working from 6 +// std::string str_rips_persistence = test_rips_persistence(6, 0); +// TODO(VR): division by zero +// std::string str_rips_persistence = test_rips_persistence(0, 0); +// TODO(VR): is result OK of : +// test_rips_persistence_in_dimension(3, 4); + diff --git a/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt.REMOVED.git-id b/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt.REMOVED.git-id new file mode 100644 index 00000000..2dd38515 --- /dev/null +++ b/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt.REMOVED.git-id @@ -0,0 +1 @@ +ce87199d425b05f51c74cbf635870bfa5abbc7a1
\ No newline at end of file diff --git a/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt b/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt new file mode 100644 index 00000000..90756ce0 --- /dev/null +++ b/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt @@ -0,0 +1,98 @@ +0 0 0 +0 1 0.01 +0 2 0.02 +0 3 0.03 +0 4 0.04 +0 5 0.05 +0 6 0.06 +0 7 0.07 +0 8 0.08 +0 9 0.09 +1 2 1 1.12 +1 3 1 1.13 +1 4 1 1.14 +1 5 1 1.15 +1 6 1 1.16 +1 7 1 1.17 +1 8 1 1.18 +1 9 1 1.19 +1 3 2 1.23 +2 3 2 1 1.23 +1 6 2 1.26 +2 6 2 1 1.26 +1 8 2 1.28 +2 8 2 1 1.28 +1 9 2 1.29 +2 9 2 1 1.29 +1 6 3 1.36 +2 6 3 1 1.36 +2 6 3 2 1.36 +3 6 3 2 1 1.36 +1 7 3 1.37 +2 7 3 1 1.37 +1 8 3 1.38 +2 8 3 1 1.38 +2 8 3 2 1.38 +3 8 3 2 1 1.38 +1 5 4 1.45 +2 5 4 1 1.45 +1 6 4 1.46 +2 6 4 1 1.46 +1 8 4 1.48 +2 8 4 1 1.48 +1 6 5 1.56 +2 6 5 1 1.56 +2 6 5 4 1.56 +3 6 5 4 1 1.56 +1 7 5 1.57 +2 7 5 1 1.57 +1 8 5 1.58 +2 8 5 1 1.58 +2 8 5 4 1.58 +3 8 5 4 1 1.58 +1 9 5 1.59 +2 9 5 1 1.59 +1 7 6 1.67 +2 7 6 1 1.67 +2 7 6 3 1.67 +3 7 6 3 1 1.67 +2 7 6 5 1.67 +3 7 6 5 1 1.67 +1 8 6 1.68 +2 8 6 1 1.68 +2 8 6 2 1.68 +3 8 6 2 1 1.68 +2 8 6 3 1.68 +3 8 6 3 1 1.68 +3 8 6 3 2 1.68 +2 8 6 4 1.68 +3 8 6 4 1 1.68 +2 8 6 5 1.68 +3 8 6 5 1 1.68 +3 8 6 5 4 1.68 +1 9 6 1.69 +2 9 6 1 1.69 +2 9 6 2 1.69 +3 9 6 2 1 1.69 +2 9 6 5 1.69 +3 9 6 5 1 1.69 +1 8 7 1.78 +2 8 7 1 1.78 +2 8 7 3 1.78 +3 8 7 3 1 1.78 +2 8 7 5 1.78 +3 8 7 5 1 1.78 +2 8 7 6 1.78 +3 8 7 6 1 1.78 +3 8 7 6 3 1.78 +3 8 7 6 5 1.78 +1 9 8 1.89 +2 9 8 1 1.89 +2 9 8 2 1.89 +3 9 8 2 1 1.89 +2 9 8 5 1.89 +3 9 8 5 1 1.89 +2 9 8 6 1.89 +3 9 8 6 1 1.89 +3 9 8 6 2 1.89 +3 9 8 6 5 1.89 diff --git a/src/Simplex_tree/example/CMakeLists.txt b/src/Simplex_tree/example/CMakeLists.txt index 6a392d3b..ac145df5 100644 --- a/src/Simplex_tree/example/CMakeLists.txt +++ b/src/Simplex_tree/example/CMakeLists.txt @@ -1,15 +1,14 @@ cmake_minimum_required(VERSION 2.6) project(GUDHISimplexTreeFromFile) -# NEED TO FIND BOOST NEEDED COMPONENTS TO LINK WITH -find_package(Boost 1.45.0 COMPONENTS system thread) -message("Boost_lib = ${Boost_LIBRARIES}") - add_executable ( simplex_tree_from_file simplex_tree_from_file.cpp ) -target_link_libraries(simplex_tree_from_file ${Boost_LIBRARIES}) - +target_link_libraries(simplex_tree_from_file ${Boost_SYSTEM_LIBRARY}) +add_test(simplex_tree_from_file_2 ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_from_file ${CMAKE_SOURCE_DIR}/data/points/Klein_bottle_complex.txt 2) +add_test(simplex_tree_from_file_3 ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_from_file ${CMAKE_SOURCE_DIR}/data/points/Klein_bottle_complex.txt 3) + add_executable ( simple_simplex_tree simple_simplex_tree.cpp ) -target_link_libraries(simple_simplex_tree ${Boost_LIBRARIES}) +target_link_libraries(simple_simplex_tree ${Boost_SYSTEM_LIBRARY}) +add_test(simple_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/simple_simplex_tree) # An example with Simplex-tree using CGAL alpha_shapes_3 if(GMPXX_FOUND AND GMP_FOUND AND CGAL_FOUND) @@ -18,5 +17,6 @@ if(GMPXX_FOUND AND GMP_FOUND AND CGAL_FOUND) message("GMP_LIBRARIES = ${GMP_LIBRARIES}") add_executable ( simplex_tree_from_alpha_shapes_3 simplex_tree_from_alpha_shapes_3.cpp ) - target_link_libraries(simplex_tree_from_alpha_shapes_3 ${Boost_LIBRARIES} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY}) + target_link_libraries(simplex_tree_from_alpha_shapes_3 ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY}) + add_test(simplex_tree_from_alpha_shapes_3 ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_from_alpha_shapes_3 ${CMAKE_SOURCE_DIR}/data/points/bunny_5000) endif() diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp index 0ed2039f..11b96a0d 100644 --- a/src/Simplex_tree/example/simple_simplex_tree.cpp +++ b/src/Simplex_tree/example/simple_simplex_tree.cpp @@ -48,12 +48,12 @@ int main (int argc, char * const argv[]) //Construct the Simplex Tree Simplex_tree<> simplexTree; - // Simplex to be inserted: - // 1 - // o - // /X\ - // o---o---o - // 2 0 3 + /* Simplex to be inserted: */ + /* 1 */ + /* o */ + /* /X\ */ + /* o---o---o */ + /* 2 0 3 */ // ++ FIRST std::cout << " * INSERT 0" << std::endl; @@ -299,4 +299,5 @@ int main (int argc, char * const argv[]) // [0.2] 2 1 // [0.2] 3 0 // [0.3] 2 1 0 + return 0; } diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 44f7a0d6..12be6e5d 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1,799 +1,838 @@ - /* 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): Clément Maria - * - * 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 GUDHI_SIMPLEX_TREE_H -#define GUDHI_SIMPLEX_TREE_H +/* 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): Clément Maria + * + * 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 SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_ +#define SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_ + +#include <gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h> +#include <gudhi/Simplex_tree/Simplex_tree_siblings.h> +#include <gudhi/Simplex_tree/Simplex_tree_iterators.h> +#include <gudhi/Simplex_tree/indexing_tag.h> -#include <algorithm> #include <boost/container/flat_map.hpp> #include <boost/iterator/transform_iterator.hpp> #include <boost/graph/adjacency_list.hpp> -#include "gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h" -#include "gudhi/Simplex_tree/Simplex_tree_siblings.h" -#include "gudhi/Simplex_tree/Simplex_tree_iterators.h" -#include "gudhi/Simplex_tree/indexing_tag.h" - -namespace Gudhi{ - -/** \defgroup simplex_tree Filtered Complexes - * - * A simplicial complex \f$\mathbf{K}\f$ - * on a set of vertices \f$V = \{1, \cdots ,|V|\}\f$ is a collection of simplices - * \f$\{\sigma\}\f$, - * \f$\sigma \subseteq V\f$ such that \f$\tau \subseteq \sigma \in \mathbf{K} \rightarrow \tau \in - * \mathbf{K}\f$. The - * dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$ is its number of elements minus \f$1\f$. - * - * A filtration of a simplicial complex is - * a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq f(\sigma)\f$ whenever - * \f$\tau \subseteq \sigma\f$. Ordering the simplices by increasing filtration values - * (breaking ties so as a simplex appears after its subsimplices of same filtration value) - * provides an indexing scheme. - * - - <DT>Implementations:</DT> - There are two implementation of complexes. The first on is the Simplex_tree data structure. - The simplex tree is an efficient and flexible - data structure for representing general (filtered) simplicial complexes. The data structure - is described in \cite boissonnatmariasimplextreealgorithmica - - The second one is the Hasse_complex. The Hasse complex is a data structure representing - explicitly all co-dimension 1 incidence relations in a complex. It is consequently faster - when accessing the boundary of a simplex, but is less compact and harder to construct from - scratch. - - - * \author Clément Maria - * \version 1.0 - * \date 2014 - * \copyright GNU General Public License v3. - * @{ - */ + +#include <algorithm> +#include <utility> +#include <vector> + +namespace Gudhi { + +/** \defgroup simplex_tree Filtered Complexes + * + * A simplicial complex \f$\mathbf{K}\f$ + * on a set of vertices \f$V = \{1, \cdots ,|V|\}\f$ is a collection of simplices + * \f$\{\sigma\}\f$, + * \f$\sigma \subseteq V\f$ such that \f$\tau \subseteq \sigma \in \mathbf{K} \rightarrow \tau \in + * \mathbf{K}\f$. The + * dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$ is its number of elements minus \f$1\f$. + * + * A filtration of a simplicial complex is + * a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq f(\sigma)\f$ whenever + * \f$\tau \subseteq \sigma\f$. Ordering the simplices by increasing filtration values + * (breaking ties so as a simplex appears after its subsimplices of same filtration value) + * provides an indexing scheme. + * + + <DT>Implementations:</DT> + There are two implementation of complexes. The first on is the Simplex_tree data structure. + The simplex tree is an efficient and flexible + data structure for representing general (filtered) simplicial complexes. The data structure + is described in \cite boissonnatmariasimplextreealgorithmica + + The second one is the Hasse_complex. The Hasse complex is a data structure representing + explicitly all co-dimension 1 incidence relations in a complex. It is consequently faster + when accessing the boundary of a simplex, but is less compact and harder to construct from + scratch. + + + * \author Clément Maria + * \version 1.0 + * \date 2014 + * \copyright GNU General Public License v3. + * @{ + */ /** - * \brief Simplex Tree data structure for representing simplicial complexes. - * - * \details Every simplex \f$[v_0, \cdots ,v_d]\f$ admits a canonical orientation - * induced by the order relation on vertices \f$ v_0 < \cdots < v_d \f$. - * - * Details may be found in \cite boissonnatmariasimplextreealgorithmica. - * - * \implements FilteredComplex - * - */ -template < typename IndexingTag = linear_indexing_tag - , typename FiltrationValue = double - , typename SimplexKey = int //must be a signed integer type - , typename VertexHandle = int //must be a signed integer type, int convertible to it + * \brief Simplex Tree data structure for representing simplicial complexes. + * + * \details Every simplex \f$[v_0, \cdots ,v_d]\f$ admits a canonical orientation + * induced by the order relation on vertices \f$ v_0 < \cdots < v_d \f$. + * + * Details may be found in \cite boissonnatmariasimplextreealgorithmica. + * + * \implements FilteredComplex + * + */ +template<typename IndexingTag = linear_indexing_tag, + typename FiltrationValue = double, typename SimplexKey = int // must be a signed integer type + , typename VertexHandle = int // must be a signed integer type, int convertible to it // , bool ContiguousVertexHandles = true //true is Vertex_handles are exactly the set [0;n) - > -class Simplex_tree -{ - -public: - typedef IndexingTag Indexing_tag; -/** \brief Type for the value of the filtration function. - * - * Must be comparable with <. */ - typedef FiltrationValue Filtration_value; -/** \brief Key associated to each simplex. - * - * Must be a signed integer type. */ - typedef SimplexKey Simplex_key; -/** \brief Type for the vertex handle. - * - * Must be a signed integer type. It admits a total order <. */ - typedef VertexHandle Vertex_handle; - - - /* Type of node in the simplex tree. */ - typedef Simplex_tree_node_explicit_storage < Simplex_tree > Node; +> +class Simplex_tree { + public: + typedef IndexingTag Indexing_tag; + /** \brief Type for the value of the filtration function. + * + * Must be comparable with <. */ + typedef FiltrationValue Filtration_value; + /** \brief Key associated to each simplex. + * + * Must be a signed integer type. */ + typedef SimplexKey Simplex_key; + /** \brief Type for the vertex handle. + * + * Must be a signed integer type. It admits a total order <. */ + typedef VertexHandle Vertex_handle; + + /* Type of node in the simplex tree. */ + typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node; /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */ - typedef typename boost::container::flat_map< Vertex_handle - , Node > Dictionary; - - - friend class Simplex_tree_node_explicit_storage < Simplex_tree < FiltrationValue - , SimplexKey - , VertexHandle > >; - friend class Simplex_tree_siblings < Simplex_tree < FiltrationValue - , SimplexKey - , VertexHandle > - , Dictionary > ; - friend class Simplex_tree_simplex_vertex_iterator < Simplex_tree < FiltrationValue - , SimplexKey - , VertexHandle > >; - friend class Simplex_tree_boundary_simplex_iterator < Simplex_tree < FiltrationValue - , SimplexKey - , VertexHandle > >; - friend class Simplex_tree_complex_simplex_iterator < Simplex_tree < FiltrationValue - , SimplexKey - , VertexHandle > >; - friend class Simplex_tree_skeleton_simplex_iterator < Simplex_tree < FiltrationValue - , SimplexKey - , VertexHandle > >; - template < class T1, class T2 > friend class Persistent_cohomology; + typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary; + + friend class Simplex_tree_node_explicit_storage< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >; + friend class Simplex_tree_siblings< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle>, Dictionary>; + friend class Simplex_tree_simplex_vertex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >; + friend class Simplex_tree_boundary_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >; + friend class Simplex_tree_complex_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >; + friend class Simplex_tree_skeleton_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >; + template<class T1, class T2> friend class Persistent_cohomology; /* \brief Set of nodes sharing a same parent in the simplex tree. */ - typedef Simplex_tree_siblings < Simplex_tree - , Dictionary > Siblings; -public: -/** \brief Handle type to a simplex contained in the simplicial complex represented - * byt he simplex tree. */ - typedef typename Dictionary::iterator Simplex_handle; -private: - typedef typename Dictionary::iterator Dictionary_it; - typedef typename Dictionary_it::value_type Dit_value_t; + typedef Simplex_tree_siblings<Simplex_tree, Dictionary> Siblings; + + public: + /** \brief Handle type to a simplex contained in the simplicial complex represented + * byt he simplex tree. */ + typedef typename Dictionary::iterator Simplex_handle; + + private: + typedef typename Dictionary::iterator Dictionary_it; + typedef typename Dictionary_it::value_type Dit_value_t; struct return_first { - Vertex_handle operator()(const Dit_value_t& p_sh) const {return p_sh.first;} + Vertex_handle operator()(const Dit_value_t& p_sh) const { + return p_sh.first; + } }; -public: -/** \name Range and iterator types - * - * The naming convention is Container_content_(iterator/range). A Container_content_range is - * essentially an object on which the methods begin() and end() can be called. They both return - * an object of type Container_content_iterator, and allow the traversal of the range - * [ begin();end() ). - * @{ */ + + public: + /** \name Range and iterator types + * + * The naming convention is Container_content_(iterator/range). A Container_content_range is + * essentially an object on which the methods begin() and end() can be called. They both return + * an object of type Container_content_iterator, and allow the traversal of the range + * [ begin();end() ). + * @{ */ /** \brief Iterator over the vertices of the simplicial complex. - * - * 'value_type' is Vertex_handle. */ - typedef boost::transform_iterator < return_first, Dictionary_it > Complex_vertex_iterator ; + * + * 'value_type' is Vertex_handle. */ + typedef boost::transform_iterator<return_first, Dictionary_it> Complex_vertex_iterator; /** \brief Range over the vertices of the simplicial complex. */ - typedef boost::iterator_range < Complex_vertex_iterator > Complex_vertex_range ; + typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range; /** \brief Iterator over the vertices of a simplex. - * - * 'value_type' is Vertex_handle. */ - typedef Simplex_tree_simplex_vertex_iterator < Simplex_tree > Simplex_vertex_iterator ; + * + * 'value_type' is Vertex_handle. */ + typedef Simplex_tree_simplex_vertex_iterator<Simplex_tree> Simplex_vertex_iterator; /** \brief Range over the vertices of a simplex. */ - typedef boost::iterator_range < Simplex_vertex_iterator > Simplex_vertex_range ; + typedef boost::iterator_range<Simplex_vertex_iterator> Simplex_vertex_range; /** \brief Iterator over the simplices of the boundary of a simplex. - * - * 'value_type' is Simplex_handle. */ - typedef Simplex_tree_boundary_simplex_iterator < Simplex_tree > Boundary_simplex_iterator ; + * + * 'value_type' is Simplex_handle. */ + typedef Simplex_tree_boundary_simplex_iterator<Simplex_tree> Boundary_simplex_iterator; /** \brief Range over the simplices of the boundary of a simplex. */ - typedef boost::iterator_range < Boundary_simplex_iterator > Boundary_simplex_range ; + typedef boost::iterator_range<Boundary_simplex_iterator> Boundary_simplex_range; /** \brief Iterator over the simplices of the simplicial complex. - * - * 'value_type' is Simplex_handle. */ - typedef Simplex_tree_complex_simplex_iterator < Simplex_tree > Complex_simplex_iterator ; + * + * 'value_type' is Simplex_handle. */ + typedef Simplex_tree_complex_simplex_iterator<Simplex_tree> Complex_simplex_iterator; /** \brief Range over the simplices of the simplicial complex. */ - typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range ; + typedef boost::iterator_range<Complex_simplex_iterator> Complex_simplex_range; /** \brief Iterator over the simplices of the skeleton of the simplicial complex, for a given - * dimension. - * - * 'value_type' is Simplex_handle. */ - typedef Simplex_tree_skeleton_simplex_iterator < Simplex_tree > Skeleton_simplex_iterator ; + * dimension. + * + * 'value_type' is Simplex_handle. */ + typedef Simplex_tree_skeleton_simplex_iterator<Simplex_tree> Skeleton_simplex_iterator; /** \brief Range over the simplices of the skeleton of the simplicial complex, for a given - * dimension. */ - typedef boost::iterator_range < Skeleton_simplex_iterator > Skeleton_simplex_range ; + * dimension. */ + typedef boost::iterator_range<Skeleton_simplex_iterator> Skeleton_simplex_range; /** \brief Iterator over the simplices of the simplicial complex, ordered by the filtration. - * - * 'value_type' is Simplex_handle. */ - typedef typename std::vector < Simplex_handle >::iterator Filtration_simplex_iterator; - /** \brief Range over the simplices of the simplicial complex, ordered by the filtration. */ - typedef boost::iterator_range < Filtration_simplex_iterator > Filtration_simplex_range ; - -/* @} */ //end name range and iterator types - - -/** \name Range and iterator methods - * @{ */ - -/** \brief Returns a range over the vertices of the simplicial complex. - * - * The order is increasing according to < on Vertex_handles.*/ -Complex_vertex_range complex_vertex_range() -{ return Complex_vertex_range( boost::make_transform_iterator(root_.members_.begin(),return_first()) - , boost::make_transform_iterator(root_.members_.end(),return_first())); } -/** \brief Returns a range over the simplices of the simplicial complex. - * - * In the Simplex_tree, the tree is traverse in a depth-first fashion. - * Consequently, simplices are ordered according to lexicographic order on the list of - * Vertex_handles of a simplex, read in increasing < order for Vertex_handles. */ - Complex_simplex_range complex_simplex_range() - { return Complex_simplex_range ( Complex_simplex_iterator(this), - Complex_simplex_iterator() ); } -/** \brief Returns a range over the simplices of the dim-skeleton of the simplicial complex. - * - * The \f$d\f$-skeleton of a simplicial complex \f$\mathbf{K}\f$ is the simplicial complex containing the - * simplices of \f$\mathbf{K}\f$ of dimension at most \f$d\f$. - * - * @param[in] dim The maximal dimension of the simplices in the skeleton. - * - * The simplices are ordered according to lexicographic order on the list of - * Vertex_handles of a simplex, read in increasing < order for Vertex_handles. */ - Skeleton_simplex_range skeleton_simplex_range(int dim) - { return Skeleton_simplex_range ( Skeleton_simplex_iterator(this,dim) - , Skeleton_simplex_iterator() ); } -/** \brief Returns a range over the simplices of the simplicial complex, - * in the order of the filtration. - * - * The filtration is a monotonic function \f$ f: \mathbf{K} \rightarrow \mathbb{R} \f$, i.e. if two simplices - * \f$\tau\f$ and \f$\sigma\f$ satisfy \f$\tau \subseteq \sigma\f$ then - * \f$f(\tau) \leq f(\sigma)\f$. - * - * The method returns simplices ordered according to increasing filtration values. Ties are - * resolved by considering inclusion relation (subsimplices appear before their cofaces). If two - * simplices have same filtration value but are not comparable w.r.t. inclusion, lexicographic - * order is used. - * - * The filtration must be valid. If the filtration has not been initialized yet, the - * method initializes it (i.e. order the simplices). */ - Filtration_simplex_range filtration_simplex_range(linear_indexing_tag) - { if(filtration_vect_.empty()) { initialize_filtration(); } - return Filtration_simplex_range ( filtration_vect_.begin() - , filtration_vect_.end()); } + * + * 'value_type' is Simplex_handle. */ + typedef typename std::vector<Simplex_handle>::iterator Filtration_simplex_iterator; + /** \brief Range over the simplices of the simplicial complex, ordered by the filtration. */ + typedef boost::iterator_range<Filtration_simplex_iterator> Filtration_simplex_range; + + /* @} */ // end name range and iterator types + /** \name Range and iterator methods + * @{ */ + + /** \brief Returns a range over the vertices of the simplicial complex. + * + * The order is increasing according to < on Vertex_handles.*/ + Complex_vertex_range complex_vertex_range() { + return Complex_vertex_range( + boost::make_transform_iterator(root_.members_.begin(), return_first()), + boost::make_transform_iterator(root_.members_.end(), return_first())); + } + + /** \brief Returns a range over the simplices of the simplicial complex. + * + * In the Simplex_tree, the tree is traverse in a depth-first fashion. + * Consequently, simplices are ordered according to lexicographic order on the list of + * Vertex_handles of a simplex, read in increasing < order for Vertex_handles. */ + Complex_simplex_range complex_simplex_range() { + return Complex_simplex_range(Complex_simplex_iterator(this), + Complex_simplex_iterator()); + } + + /** \brief Returns a range over the simplices of the dim-skeleton of the simplicial complex. + * + * The \f$d\f$-skeleton of a simplicial complex \f$\mathbf{K}\f$ is the simplicial complex containing the + * simplices of \f$\mathbf{K}\f$ of dimension at most \f$d\f$. + * + * @param[in] dim The maximal dimension of the simplices in the skeleton. + * + * The simplices are ordered according to lexicographic order on the list of + * Vertex_handles of a simplex, read in increasing < order for Vertex_handles. */ + Skeleton_simplex_range skeleton_simplex_range(int dim) { + return Skeleton_simplex_range(Skeleton_simplex_iterator(this, dim), + Skeleton_simplex_iterator()); + } + + /** \brief Returns a range over the simplices of the simplicial complex, + * in the order of the filtration. + * + * The filtration is a monotonic function \f$ f: \mathbf{K} \rightarrow \mathbb{R} \f$, i.e. if two simplices + * \f$\tau\f$ and \f$\sigma\f$ satisfy \f$\tau \subseteq \sigma\f$ then + * \f$f(\tau) \leq f(\sigma)\f$. + * + * The method returns simplices ordered according to increasing filtration values. Ties are + * resolved by considering inclusion relation (subsimplices appear before their cofaces). If two + * simplices have same filtration value but are not comparable w.r.t. inclusion, lexicographic + * order is used. + * + * The filtration must be valid. If the filtration has not been initialized yet, the + * method initializes it (i.e. order the simplices). */ + Filtration_simplex_range filtration_simplex_range(linear_indexing_tag) { + if (filtration_vect_.empty()) { + initialize_filtration(); + } + return Filtration_simplex_range(filtration_vect_.begin(), + filtration_vect_.end()); + } Filtration_simplex_range filtration_simplex_range() { return filtration_simplex_range(Indexing_tag()); } -/** \brief Returns a range over the vertices of a simplex. - * - * The order in which the vertices are visited is the decreasing order for < on Vertex_handles, - * which is consequenlty - * equal to \f$(-1)^{\text{dim} \sigma}\f$ the canonical orientation on the simplex. - */ - Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) - { return Simplex_vertex_range ( Simplex_vertex_iterator(this,sh), - Simplex_vertex_iterator(this));} - -/** \brief Returns a range over the simplices of the boundary of a simplex. - * - * The boundary of a simplex is the set of codimension \f$1\f$ subsimplices of the simplex. - * If the simplex is \f$[v_0, \cdots ,v_d]\f$, with canonical orientation - * induced by \f$ v_0 < \cdots < v_d \f$, the iterator enumerates the - * simplices of the boundary in the order: - * \f$[v_0,\cdots,\widehat{v_i},\cdots,v_d]\f$ for \f$i\f$ from \f$0\f$ to \f$d\f$, - * where \f$\widehat{v_i}\f$ means that the vertex \f$v_i\f$ is omitted. - * - * We note that the alternate sum of the simplices given by the iterator - * gives \f$(-1)^{\text{dim} \sigma}\f$ the chains corresponding to the boundary - * of the simplex. - * - * @param[in] sh Simplex for which the boundary is computed. */ - Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) - { return Boundary_simplex_range ( Boundary_simplex_iterator(this,sh), - Boundary_simplex_iterator(this) ); } - -/** @} */ //end range and iterator methods - - -/** \name Constructor/Destructor - * @{ */ - -/** \brief Constructs an empty simplex tree. */ - Simplex_tree () - : null_vertex_(-1) - , threshold_(0) - , num_simplices_(0) - , root_(NULL,null_vertex_) - , filtration_vect_() - , dimension_(-1) {} - -/** \brief Destructor; deallocates the whole tree structure. */ -~Simplex_tree() -{ - for(auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) - { - if(has_children(sh)) { rec_delete(sh->second.children()); } - } -} -/** @} */ // end constructor/destructor -private: -/** Recursive deletion. */ -void rec_delete(Siblings * sib) -{ - for(auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) - { if(has_children(sh)) { rec_delete(sh->second.children()); } } - delete sib; -} + /** \brief Returns a range over the vertices of a simplex. + * + * The order in which the vertices are visited is the decreasing order for < on Vertex_handles, + * which is consequenlty + * equal to \f$(-1)^{\text{dim} \sigma}\f$ the canonical orientation on the simplex. + */ + Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) { + return Simplex_vertex_range(Simplex_vertex_iterator(this, sh), + Simplex_vertex_iterator(this)); + } + /** \brief Returns a range over the simplices of the boundary of a simplex. + * + * The boundary of a simplex is the set of codimension \f$1\f$ subsimplices of the simplex. + * If the simplex is \f$[v_0, \cdots ,v_d]\f$, with canonical orientation + * induced by \f$ v_0 < \cdots < v_d \f$, the iterator enumerates the + * simplices of the boundary in the order: + * \f$[v_0,\cdots,\widehat{v_i},\cdots,v_d]\f$ for \f$i\f$ from \f$0\f$ to \f$d\f$, + * where \f$\widehat{v_i}\f$ means that the vertex \f$v_i\f$ is omitted. + * + * We note that the alternate sum of the simplices given by the iterator + * gives \f$(-1)^{\text{dim} \sigma}\f$ the chains corresponding to the boundary + * of the simplex. + * + * @param[in] sh Simplex for which the boundary is computed. */ + Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { + return Boundary_simplex_range(Boundary_simplex_iterator(this, sh), + Boundary_simplex_iterator(this)); + } -public: -/** \brief Returns the key associated to a simplex. - * - * The filtration must be initialized. */ -Simplex_key key ( Simplex_handle sh ) { return sh->second.key(); } -/** \brief Returns the simplex associated to a key. - * - * The filtration must be initialized. */ -Simplex_handle simplex ( Simplex_key key ) { return filtration_vect_[key]; } -/** \brief Returns the filtration value of a simplex. - * - * Called on the null_simplex, returns INFINITY. */ -Filtration_value filtration(Simplex_handle sh) -{ - if(sh != null_simplex()) { return sh->second.filtration(); } - else { return INFINITY; }//filtration(); } -} -/** \brief Returns an upper bound of the filtration values of the simplices. */ -Filtration_value filtration() -{ return threshold_; } -/** \brief Returns a Simplex_handle different from all Simplex_handles - * associated to the simplices in the simplicial complex. - * - * One can call filtration(null_simplex()). */ -Simplex_handle null_simplex() { return Dictionary_it(NULL); } -/** \brief Returns a key different for all keys associated to the - * simplices of the simplicial complex. */ -Simplex_key null_key() { return -1; } -/** \brief Returns a Vertex_handle different from all Vertex_handles associated - * to the vertices of the simplicial complex. */ -Vertex_handle null_vertex() { return null_vertex_; } -/** \brief Returns the number of vertices in the complex. */ -size_t num_vertices() { return root_.members_.size(); } -/** \brief Returns the number of simplices in the complex. - * - * Does not count the empty simplex. */ -size_t num_simplices() { return num_simplices_; } - -/** \brief Returns the dimension of a simplex. - * - * Must be different from null_simplex().*/ -int dimension(Simplex_handle sh) -{ Siblings * curr_sib = self_siblings(sh); - int dim = 0; - while(curr_sib != NULL) { ++dim; curr_sib = curr_sib->oncles(); } - return dim-1; -} -/** \brief Returns an upper bound on the dimension of the simplicial complex. */ -int dimension() -{ return dimension_; } - -/** \brief Returns true iff the node in the simplex tree pointed by - * sh has children.*/ -bool has_children(Simplex_handle sh) -{ return (sh->second.children()->parent() == sh->first); } - -/** \brief Given a range of Vertex_handles, returns the Simplex_handle - * of the simplex in the simplicial complex containing the corresponding - * vertices. Return null_simplex() if the simplex is not in the complex. - * - * The type RandomAccessVertexRange must be a range for which .begin() and - * .end() return random access iterators, with <CODE>value_type</CODE> - * <CODE>Vertex_handle</CODE>. - */ -template <class RandomAccessVertexRange > -Simplex_handle find(RandomAccessVertexRange & s) -{ - if(s.begin() == s.end()) std::cerr << "Empty simplex \n"; - - sort(s.begin(),s.end()); - - Siblings * tmp_sib = &root_; - Dictionary_it tmp_dit; - Vertex_handle last = s[s.size()-1]; - for(auto v : s) { - tmp_dit = tmp_sib->members_.find(v); - if(tmp_dit == tmp_sib->members_.end()) { return null_simplex(); } - if( !has_children(tmp_dit) && v != last) { return null_simplex(); } - tmp_sib = tmp_dit->second.children(); - } - return tmp_dit; -} + /** @} */ // end range and iterator methods + /** \name Constructor/Destructor + * @{ */ + + /** \brief Constructs an empty simplex tree. */ + Simplex_tree() + : null_vertex_(-1), + threshold_(0), + num_simplices_(0), + root_(NULL, null_vertex_), + filtration_vect_(), + dimension_(-1) { + } + + /** \brief Destructor; deallocates the whole tree structure. */ + ~Simplex_tree() { + for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) { + if (has_children(sh)) { + rec_delete(sh->second.children()); + } + } + } + /** @} */ // end constructor/destructor + private: + /** Recursive deletion. */ + void rec_delete(Siblings * sib) { + for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) { + if (has_children(sh)) { + rec_delete(sh->second.children()); + } + } + delete sib; + } + + public: + /** \brief Returns the key associated to a simplex. + * + * The filtration must be initialized. */ + Simplex_key key(Simplex_handle sh) { + return sh->second.key(); + } + /** \brief Returns the simplex associated to a key. + * + * The filtration must be initialized. */ + Simplex_handle simplex(Simplex_key key) { + return filtration_vect_[key]; + } + /** \brief Returns the filtration value of a simplex. + * + * Called on the null_simplex, returns INFINITY. */ + Filtration_value filtration(Simplex_handle sh) { + if (sh != null_simplex()) { + return sh->second.filtration(); + } else { + return INFINITY; + } // filtration(); } + } + /** \brief Returns an upper bound of the filtration values of the simplices. */ + Filtration_value filtration() { + return threshold_; + } + /** \brief Returns a Simplex_handle different from all Simplex_handles + * associated to the simplices in the simplicial complex. + * + * One can call filtration(null_simplex()). */ + Simplex_handle null_simplex() { + return Dictionary_it(NULL); + } + /** \brief Returns a key different for all keys associated to the + * simplices of the simplicial complex. */ + Simplex_key null_key() { + return -1; + } + /** \brief Returns a Vertex_handle different from all Vertex_handles associated + * to the vertices of the simplicial complex. */ + Vertex_handle null_vertex() { + return null_vertex_; + } + /** \brief Returns the number of vertices in the complex. */ + size_t num_vertices() { + return root_.members_.size(); + } + /** \brief Returns the number of simplices in the complex. + * + * Does not count the empty simplex. */ + const unsigned int& num_simplices() const { + return num_simplices_; + } + + /** \brief Returns the dimension of a simplex. + * + * Must be different from null_simplex().*/ + int dimension(Simplex_handle sh) { + Siblings * curr_sib = self_siblings(sh); + int dim = 0; + while (curr_sib != NULL) { + ++dim; + curr_sib = curr_sib->oncles(); + } + return dim - 1; + } + /** \brief Returns an upper bound on the dimension of the simplicial complex. */ + int dimension() { + return dimension_; + } + + /** \brief Returns true iff the node in the simplex tree pointed by + * sh has children.*/ + bool has_children(Simplex_handle sh) { + return (sh->second.children()->parent() == sh->first); + } -/** \brief Returns the Simplex_handle corresponding to the 0-simplex - * representing the vertex with Vertex_handle v. */ -Simplex_handle find_vertex(Vertex_handle v) -{ return root_.members_.begin()+v; } + /** \brief Given a range of Vertex_handles, returns the Simplex_handle + * of the simplex in the simplicial complex containing the corresponding + * vertices. Return null_simplex() if the simplex is not in the complex. + * + * The type RandomAccessVertexRange must be a range for which .begin() and + * .end() return random access iterators, with <CODE>value_type</CODE> + * <CODE>Vertex_handle</CODE>. + */ + template<class RandomAccessVertexRange> + Simplex_handle find(const RandomAccessVertexRange & s) { + if (s.begin() == s.end()) + std::cerr << "Empty simplex \n"; + + sort(s.begin(), s.end()); + + Siblings * tmp_sib = &root_; + Dictionary_it tmp_dit; + Vertex_handle last = s[s.size() - 1]; + for (auto v : s) { + tmp_dit = tmp_sib->members_.find(v); + if (tmp_dit == tmp_sib->members_.end()) { + return null_simplex(); + } + if (!has_children(tmp_dit) && v != last) { + return null_simplex(); + } + tmp_sib = tmp_dit->second.children(); + } + return tmp_dit; + } + + /** \brief Returns the Simplex_handle corresponding to the 0-simplex + * representing the vertex with Vertex_handle v. */ + Simplex_handle find_vertex(Vertex_handle v) { + return root_.members_.begin() + v; + } //{ return root_.members_.find(v); } + /** \brief Insert a simplex, represented by a range of Vertex_handles, in the simplicial complex. + * + * @param[in] simplex range of Vertex_handles, representing the vertices of the new simplex + * @param[in] filtration the filtration value assigned to the new simplex. + * The return type is a pair. If the new simplex is inserted successfully (i.e. it was not in the + * simplicial complex yet) the bool is set to true and the Simplex_handle is the handle assigned + * to the new simplex. + * If the insertion fails (the simplex is already there), the bool is set to false. If the insertion + * fails and the simplex already in the complex has a filtration value strictly bigger than 'filtration', + * we assign this simplex with the new value 'filtration', and set the Simplex_handle filed of the + * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to + * null_simplex. + * + * All subsimplices do not necessary need to be already in the simplex tree to proceed to an + * insertion. However, the property of being a simplicial complex will be violated. This allows + * us to insert a stream of simplices contained in a simplicial complex without considering any + * order on them. + * + * The filtration value + * assigned to the new simplex must preserve the monotonicity of the filtration. + * + * The type RandomAccessVertexRange must be a range for which .begin() and + * .end() return random access iterators, with 'value_type' Vertex_handle. */ + template<class RandomAccessVertexRange> + std::pair<Simplex_handle, bool> insert(RandomAccessVertexRange & simplex, + Filtration_value filtration) { + if (simplex.empty()) { + return std::pair<Simplex_handle, bool>(null_simplex(), true); + } + + sort(simplex.begin(), simplex.end()); // must be sorted in increasing order -/** \brief Insert a simplex, represented by a range of Vertex_handles, in the simplicial complex. - * - * @param[in] simplex range of Vertex_handles, representing the vertices of the new simplex - * @param[in] filtration the filtration value assigned to the new simplex. - * The return type is a pair. If the new simplex is inserted successfully (i.e. it was not in the - * simplicial complex yet) the bool is set to true and the Simplex_handle is the handle assigned - * to the new simplex. - * If the insertion fails (the simplex is already there), the bool is set to false. If the insertion - * fails and the simplex already in the complex has a filtration value strictly bigger than 'filtration', - * we assign this simplex with the new value 'filtration', and set the Simplex_handle filed of the - * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to - * null_simplex. - * - * All subsimplices do not necessary need to be already in the simplex tree to proceed to an - * insertion. However, the property of being a simplicial complex will be violated. This allows - * us to insert a stream of simplices contained in a simplicial complex without considering any - * order on them. - * - * The filtration value - * assigned to the new simplex must preserve the monotonicity of the filtration. - * - * The type RandomAccessVertexRange must be a range for which .begin() and - * .end() return random access iterators, with 'value_type' Vertex_handle. */ -template <class RandomAccessVertexRange > -std::pair< Simplex_handle, bool > -insert ( RandomAccessVertexRange & simplex - , Filtration_value filtration ) -{ - if(simplex.empty()) { return std::pair< Simplex_handle, bool >(null_simplex(),true); } - - sort(simplex.begin(),simplex.end()); //must be sorted in increasing order - - Siblings * curr_sib = &root_; - std::pair< Simplex_handle, bool > res_insert; - typename RandomAccessVertexRange::iterator vi; - for( vi = simplex.begin(); vi != simplex.end()-1; ++vi ) - { - res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,filtration)); - if(!(has_children(res_insert.first))) - { res_insert.first->second.assign_children( new Siblings(curr_sib, *vi)); } - curr_sib = res_insert.first->second.children(); - } - res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,filtration)); - if(!res_insert.second) //if already in the complex - { - if(res_insert.first->second.filtration() > filtration) //if filtration value modified - { + Siblings * curr_sib = &root_; + std::pair<Simplex_handle, bool> res_insert; + typename RandomAccessVertexRange::iterator vi; + for (vi = simplex.begin(); vi != simplex.end() - 1; ++vi) { + res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration)); + if (!(has_children(res_insert.first))) { + res_insert.first->second.assign_children(new Siblings(curr_sib, *vi)); + } + curr_sib = res_insert.first->second.children(); + } + res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration)); + if (!res_insert.second) { // if already in the complex + if (res_insert.first->second.filtration() > filtration) { // if filtration value modified res_insert.first->second.assign_filtration(filtration); return res_insert; } - return std::pair< Simplex_handle, bool > (null_simplex(),false);// if filtration value unchanged + return std::pair<Simplex_handle, bool>(null_simplex(), false); // if filtration value unchanged } - //otherwise the insertion has succeeded - return res_insert; -} + // otherwise the insertion has succeeded + return res_insert; + } -/** \brief Assign a value 'key' to the key of the simplex - * represented by the Simplex_handle 'sh'. */ -void assign_key(Simplex_handle sh, Simplex_key key) -{ sh->second.assign_key(key);} + /** \brief Assign a value 'key' to the key of the simplex + * represented by the Simplex_handle 'sh'. */ + void assign_key(Simplex_handle sh, Simplex_key key) { + sh->second.assign_key(key); + } -public: -/** Returns the two Simplex_handle corresponding to the endpoints of - * and edge. sh must point to a 1-dimensional simplex. This is an - * optimized version of the boundary computation. */ -std::pair<Simplex_handle,Simplex_handle> endpoints(Simplex_handle sh) -{ return std::pair<Simplex_handle,Simplex_handle>(root_.members_.find(sh->first) - , root_.members_.find(self_siblings(sh)->parent()) ); } + public: + /** Returns the two Simplex_handle corresponding to the endpoints of + * and edge. sh must point to a 1-dimensional simplex. This is an + * optimized version of the boundary computation. */ + std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) { + return std::pair<Simplex_handle, Simplex_handle>( + root_.members_.find(sh->first), + root_.members_.find(self_siblings(sh)->parent())); + } -/** Returns the Siblings containing a simplex.*/ -Siblings * self_siblings(Simplex_handle sh) -{ if(sh->second.children()->parent() == sh->first) return sh->second.children()->oncles(); - else return sh->second.children(); } + /** Returns the Siblings containing a simplex.*/ + Siblings * self_siblings(Simplex_handle sh) { + if (sh->second.children()->parent() == sh->first) + return sh->second.children()->oncles(); + else + return sh->second.children(); + } // void display_simplex(Simplex_handle sh) // { -// std::cout << " " << "[" << filtration(sh) << "] "; -// for( auto vertex : simplex_vertex_range(sh) ) -// { std::cout << vertex << " "; } +// std::cout << " " << "[" << filtration(sh) << "] "; +// for( auto vertex : simplex_vertex_range(sh) ) +// { std::cout << vertex << " "; } // } - // void print(Simplex_handle sh, std::ostream& os = std::cout) - // { for(auto v : simplex_vertex_range(sh)) {os << v << " ";} - // os << std::endl; } - -public: -/** Returns a pointer to the root nodes of the simplex tree. */ -Siblings * root() { return &root_; } - - - - -public: -/** Set an upper bound for the filtration values. */ -void set_filtration(Filtration_value fil) -{ threshold_ = fil; } -/** Set a number of simplices for the simplicial complex. */ -void set_num_simplices(size_t num_simplices) -{ num_simplices_ = num_simplices; } -/** Set a dimension for the simplicial complex. */ -void set_dimension(int dimension) -{ dimension_ = dimension; } - -public: -/** \brief Initializes the filtrations, i.e. sort the - * simplices according to their order in the filtration and initializes all Simplex_keys. - * - * After calling this method, filtration_simplex_range() becomes valid, and each simplex is - * assigned a Simplex_key corresponding to its order in the filtration (from 0 to m-1 for a - * simplicial complex with m simplices). - * - * The use of a depth-first traversal of the simplex tree, provided by - * complex_simplex_range(), combined with - * a stable sort is meant to optimize the order of simplices with same - * filtration value. The heuristic consists in inserting the cofaces of a - * simplex as soon as possible. - * - * Will be automatically called when calling filtration_simplex_range() - * if the filtration has never been initialized yet. */ -void initialize_filtration() -{ - filtration_vect_.clear(); - filtration_vect_.reserve(num_simplices()); - for(auto cpx_it = complex_simplex_range().begin(); - cpx_it != complex_simplex_range().end(); ++cpx_it) { filtration_vect_.push_back(*cpx_it); } - -//the stable sort ensures the ordering heuristic - std::stable_sort(filtration_vect_.begin(),filtration_vect_.end(),is_before_in_filtration(this)); -} + // void print(Simplex_handle sh, std::ostream& os = std::cout) + // { for(auto v : simplex_vertex_range(sh)) {os << v << " ";} + // os << std::endl; } -private: -/** \brief Returns true iff the list of vertices of sh1 - * is smaller than the list of vertices of sh2 w.r.t. - * lexicographic order on the lists read in reverse. - * - * It defines a StrictWeakOrdering on simplices. The Simplex_vertex_iterators - * must traverse the Vertex_handle in decreasing order. Reverse lexicographic order satisfy - * the property that a subsimplex of a simplex is always strictly smaller with this order. */ -bool reverse_lexicographic_order(Simplex_handle sh1, Simplex_handle sh2) -{ - Simplex_vertex_range rg1 = simplex_vertex_range(sh1); - Simplex_vertex_range rg2 = simplex_vertex_range(sh2); - Simplex_vertex_iterator it1 = rg1.begin(); - Simplex_vertex_iterator it2 = rg2.begin(); - while(it1 != rg1.end() && it2 != rg2.end()) - { - if(*it1 == *it2) {++it1; ++it2;} - else { return *it1 < *it2; } - } - return ( (it1 == rg1.end()) && (it2 != rg2.end()) ); -} -/** \brief StrictWeakOrdering, for the simplices, defined by the filtration. - * - * It corresponds to the partial order - * induced by the filtration values, with ties resolved using reverse lexicographic order. - * Reverse lexicographic order has the property to always consider the subsimplex of a simplex - * to be smaller. The filtration function must be monotonic. */ -struct is_before_in_filtration { - is_before_in_filtration(Simplex_tree * st) : st_(st) {} - - bool operator()( const Simplex_handle sh1, - const Simplex_handle sh2 ) const - { - if(st_->filtration(sh1) != st_->filtration(sh2)) - { return st_->filtration(sh1) < st_->filtration(sh2); } - - return st_->reverse_lexicographic_order(sh1,sh2); //is sh1 a proper subface of sh2 - } - - Simplex_tree * st_; -}; + public: + /** Returns a pointer to the root nodes of the simplex tree. */ + Siblings * root() { + return &root_; + } -public: -/** \brief Inserts a 1-skeleton in an empty Simplex_tree. - * - * The Simplex_tree must contain no simplex when the method is - * called. - * - * Inserts all vertices and edges given by a OneSkeletonGraph. - * OneSkeletonGraph must be a model of boost::AdjacencyGraph, - * boost::EdgeListGraph and boost::PropertyGraph. - * - * The vertex filtration value is accessible through the property tag - * vertex_filtration_t. - * The edge filtration value is accessible through the property tag - * edge_filtration_t. - * - * boost::graph_traits<OneSkeletonGraph>::vertex_descriptor - * must be Vertex_handle. - * boost::graph_traits<OneSkeletonGraph>::directed_category - * must be undirected_tag. */ - template< class OneSkeletonGraph > - void insert_graph( OneSkeletonGraph & skel_graph ) - { - assert(num_simplices() == 0); //the simplex tree must be empty - - if(boost::num_vertices(skel_graph) == 0) { return; } - if(num_edges(skel_graph) == 0) { dimension_ = 0; } - else { dimension_ = 1; } - - num_simplices_ = boost::num_vertices(skel_graph) + boost::num_edges(skel_graph); - root_.members_.reserve(boost::num_vertices(skel_graph)); - - typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it, v_it_end; - for(tie(v_it,v_it_end) = boost::vertices(skel_graph); - v_it != v_it_end; ++v_it) - { - root_.members_.emplace_hint( root_.members_.end() - , *v_it - , Node(&root_ - ,boost::get( vertex_filtration_t() - , skel_graph - , *v_it) ) ); - } - typename boost::graph_traits<OneSkeletonGraph>::edge_iterator e_it, e_it_end; - for(tie(e_it,e_it_end) = boost::edges(skel_graph); - e_it != e_it_end; ++e_it) - { - auto u = source( *e_it, skel_graph ); auto v = target( *e_it, skel_graph ); - if( u < v ) // count edges only once { std::swap(u,v); } // u < v - { - auto sh = find_vertex(u); - if(! has_children(sh) ) { sh->second.assign_children(new Siblings(&root_,sh->first)); } - - sh->second.children()->members().emplace ( v - , Node( sh->second.children() - , boost::get( edge_filtration_t() - , skel_graph - , *e_it) )); + public: + /** Set an upper bound for the filtration values. */ + void set_filtration(Filtration_value fil) { + threshold_ = fil; + } + /** Set a number of simplices for the simplicial complex. */ + void set_num_simplices(const unsigned int& num_simplices) { + num_simplices_ = num_simplices; + } + /** Set a dimension for the simplicial complex. */ + void set_dimension(int dimension) { + dimension_ = dimension; + } + + public: + /** \brief Initializes the filtrations, i.e. sort the + * simplices according to their order in the filtration and initializes all Simplex_keys. + * + * After calling this method, filtration_simplex_range() becomes valid, and each simplex is + * assigned a Simplex_key corresponding to its order in the filtration (from 0 to m-1 for a + * simplicial complex with m simplices). + * + * The use of a depth-first traversal of the simplex tree, provided by + * complex_simplex_range(), combined with + * a stable sort is meant to optimize the order of simplices with same + * filtration value. The heuristic consists in inserting the cofaces of a + * simplex as soon as possible. + * + * Will be automatically called when calling filtration_simplex_range() + * if the filtration has never been initialized yet. */ + void initialize_filtration() { + filtration_vect_.clear(); + filtration_vect_.reserve(num_simplices()); + for (auto cpx_it = complex_simplex_range().begin(); + cpx_it != complex_simplex_range().end(); ++cpx_it) { + filtration_vect_.push_back(*cpx_it); } + + // the stable sort ensures the ordering heuristic + std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(), + is_before_in_filtration(this)); } -} -/** \brief Expands the Simplex_tree containing only its one skeleton - * until dimension max_dim. - * - * The expanded simplicial complex until dimension \f$d\f$ - * attached to a graph \f$G\f$ is the maximal simplicial complex of - * dimension at most \f$d\f$ admitting the graph \f$G\f$ as \f$1\f$-skeleton. - * The filtration value assigned to a simplex is the maximal filtration - * value of one of its edges. - * - * The Simplex_tree must contain no simplex of dimension bigger than - * 1 when calling the method. */ -void expansion(int max_dim) -{ - dimension_ = max_dim; - for(Dictionary_it root_it = root_.members_.begin(); - root_it != root_.members_.end(); ++root_it) - { - if(has_children(root_it)) - { siblings_expansion(root_it->second.children(), max_dim-1); } - } - dimension_ = max_dim - dimension_; -} -private: -/** \brief Recursive expansion of the simplex tree.*/ -void siblings_expansion ( Siblings * siblings, //must contain elements - int k) -{ - if(dimension_ > k) { dimension_ = k; } - if(k == 0) return; - Dictionary_it next = siblings->members().begin(); ++next; - - static std::vector< std::pair<Vertex_handle , Node> > inter; // static, not thread-safe. - for(Dictionary_it s_h = siblings->members().begin(); - s_h != siblings->members().end(); ++s_h,++next) - { - Simplex_handle root_sh = find_vertex(s_h->first); - if(has_children(root_sh)) - { - intersection(inter, //output intersection - next, //begin - siblings->members().end(),//end - root_sh->second.children()->members().begin(), - root_sh->second.children()->members().end(), - s_h->second.filtration()); - if(inter.size() != 0) - { this->num_simplices_ += inter.size(); - Siblings * new_sib = new Siblings(siblings, //oncles - s_h->first, //parent - inter); //boost::container::ordered_unique_range_t - inter.clear(); - s_h->second.assign_children(new_sib); - siblings_expansion(new_sib,k-1); + + private: + /** \brief Returns true iff the list of vertices of sh1 + * is smaller than the list of vertices of sh2 w.r.t. + * lexicographic order on the lists read in reverse. + * + * It defines a StrictWeakOrdering on simplices. The Simplex_vertex_iterators + * must traverse the Vertex_handle in decreasing order. Reverse lexicographic order satisfy + * the property that a subsimplex of a simplex is always strictly smaller with this order. */ + bool reverse_lexicographic_order(Simplex_handle sh1, Simplex_handle sh2) { + Simplex_vertex_range rg1 = simplex_vertex_range(sh1); + Simplex_vertex_range rg2 = simplex_vertex_range(sh2); + Simplex_vertex_iterator it1 = rg1.begin(); + Simplex_vertex_iterator it2 = rg2.begin(); + while (it1 != rg1.end() && it2 != rg2.end()) { + if (*it1 == *it2) { + ++it1; + ++it2; + } else { + return *it1 < *it2; } - else { s_h->second.assign_children(siblings); //ensure the children property - inter.clear();} } + return ((it1 == rg1.end()) && (it2 != rg2.end())); } -} -/** \brief Intersects Dictionary 1 [begin1;end1) with Dictionary 2 [begin2,end2) - * and assigns the maximal possible Filtration_value to the Nodes. */ -void intersection ( std::vector< std::pair< Vertex_handle, Node > > & intersection - , Dictionary_it begin1 - , Dictionary_it end1 - , Dictionary_it begin2 - , Dictionary_it end2 - , Filtration_value filtration ) -{ - if(begin1 == end1 || begin2 == end2) return;// 0; - while( true ) - { - if( begin1->first == begin2->first ) - { - intersection.push_back(std::pair< Vertex_handle, Node >( begin1->first, - Node( NULL - , maximum( begin1->second.filtration() - , begin2->second.filtration() - , filtration ) - ) - ) - ); - ++begin1; ++begin2; - if( begin1 == end1 || begin2 == end2 ) return; + /** \brief StrictWeakOrdering, for the simplices, defined by the filtration. + * + * It corresponds to the partial order + * induced by the filtration values, with ties resolved using reverse lexicographic order. + * Reverse lexicographic order has the property to always consider the subsimplex of a simplex + * to be smaller. The filtration function must be monotonic. */ + struct is_before_in_filtration { + explicit is_before_in_filtration(Simplex_tree * st) + : st_(st) { } - else - { - if( begin1->first < begin2->first ) - { ++begin1; if(begin1 == end1) return; } - else { ++begin2; if(begin2 == end2) return;} + + bool operator()(const Simplex_handle sh1, const Simplex_handle sh2) const { + if (st_->filtration(sh1) != st_->filtration(sh2)) { + return st_->filtration(sh1) < st_->filtration(sh2); + } + + return st_->reverse_lexicographic_order(sh1, sh2); // is sh1 a proper subface of sh2 + } + + Simplex_tree * st_; + }; + + public: + /** \brief Inserts a 1-skeleton in an empty Simplex_tree. + * + * The Simplex_tree must contain no simplex when the method is + * called. + * + * Inserts all vertices and edges given by a OneSkeletonGraph. + * OneSkeletonGraph must be a model of boost::AdjacencyGraph, + * boost::EdgeListGraph and boost::PropertyGraph. + * + * The vertex filtration value is accessible through the property tag + * vertex_filtration_t. + * The edge filtration value is accessible through the property tag + * edge_filtration_t. + * + * boost::graph_traits<OneSkeletonGraph>::vertex_descriptor + * must be Vertex_handle. + * boost::graph_traits<OneSkeletonGraph>::directed_category + * must be undirected_tag. */ + template<class OneSkeletonGraph> + void insert_graph(const OneSkeletonGraph& skel_graph) { + assert(num_simplices() == 0); // the simplex tree must be empty + + if (boost::num_vertices(skel_graph) == 0) { + return; + } + if (num_edges(skel_graph) == 0) { + dimension_ = 0; + } else { + dimension_ = 1; + } + + num_simplices_ = boost::num_vertices(skel_graph) + + boost::num_edges(skel_graph); + root_.members_.reserve(boost::num_vertices(skel_graph)); + + typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it, + v_it_end; + for (tie(v_it, v_it_end) = boost::vertices(skel_graph); v_it != v_it_end; + ++v_it) { + root_.members_.emplace_hint( + root_.members_.end(), *v_it, + Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it))); + } + typename boost::graph_traits<OneSkeletonGraph>::edge_iterator e_it, + e_it_end; + for (tie(e_it, e_it_end) = boost::edges(skel_graph); e_it != e_it_end; + ++e_it) { + auto u = source(*e_it, skel_graph); + auto v = target(*e_it, skel_graph); + if (u < v) { // count edges only once { std::swap(u,v); } // u < v + auto sh = find_vertex(u); + if (!has_children(sh)) { + sh->second.assign_children(new Siblings(&root_, sh->first)); + } + + sh->second.children()->members().emplace( + v, + Node(sh->second.children(), + boost::get(edge_filtration_t(), skel_graph, *e_it))); + } } } -} -/** Maximum over 3 values.*/ -Filtration_value maximum( Filtration_value a, - Filtration_value b, - Filtration_value c ) -{ Filtration_value max = ( a < b ) ? b : a; - return ( ( max < c ) ? c : max ); } - -public: -/** \brief Write the hasse diagram of the simplicial complex in os. - * - * Each row in the file correspond to a simplex. A line is written: - * dim idx_1 ... idx_k fil where dim is the dimension of the simplex, - * idx_1 ... idx_k are the row index (starting from 0) of the simplices of the boundary - * of the simplex, and fil is its filtration value. */ -void print_hasse(std::ostream & os) -{ - os << num_simplices() << " " << std::endl; - for(auto sh : filtration_simplex_range(Indexing_tag())) - { - os << dimension(sh) << " "; - for(auto b_sh : boundary_simplex_range(sh)) { os << key(b_sh) << " "; } - os << filtration(sh) << " \n"; + /** \brief Expands the Simplex_tree containing only its one skeleton + * until dimension max_dim. + * + * The expanded simplicial complex until dimension \f$d\f$ + * attached to a graph \f$G\f$ is the maximal simplicial complex of + * dimension at most \f$d\f$ admitting the graph \f$G\f$ as \f$1\f$-skeleton. + * The filtration value assigned to a simplex is the maximal filtration + * value of one of its edges. + * + * The Simplex_tree must contain no simplex of dimension bigger than + * 1 when calling the method. */ + void expansion(int max_dim) { + dimension_ = max_dim; + for (Dictionary_it root_it = root_.members_.begin(); + root_it != root_.members_.end(); ++root_it) { + if (has_children(root_it)) { + siblings_expansion(root_it->second.children(), max_dim - 1); + } + } + dimension_ = max_dim - dimension_; } -} -private: + private: + /** \brief Recursive expansion of the simplex tree.*/ + void siblings_expansion(Siblings * siblings, // must contain elements + int k) { + if (dimension_ > k) { + dimension_ = k; + } + if (k == 0) + return; + Dictionary_it next = siblings->members().begin(); + ++next; + + static std::vector<std::pair<Vertex_handle, Node> > inter; // static, not thread-safe. + for (Dictionary_it s_h = siblings->members().begin(); + s_h != siblings->members().end(); ++s_h, ++next) { + Simplex_handle root_sh = find_vertex(s_h->first); + if (has_children(root_sh)) { + intersection( + inter, // output intersection + next, // begin + siblings->members().end(), // end + root_sh->second.children()->members().begin(), + root_sh->second.children()->members().end(), + s_h->second.filtration()); + if (inter.size() != 0) { + this->num_simplices_ += inter.size(); + Siblings * new_sib = new Siblings(siblings, // oncles + s_h->first, // parent + inter); // boost::container::ordered_unique_range_t + inter.clear(); + s_h->second.assign_children(new_sib); + siblings_expansion(new_sib, k - 1); + } else { + s_h->second.assign_children(siblings); // ensure the children property + inter.clear(); + } + } + } + } + /** \brief Intersects Dictionary 1 [begin1;end1) with Dictionary 2 [begin2,end2) + * and assigns the maximal possible Filtration_value to the Nodes. */ + void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection, + Dictionary_it begin1, Dictionary_it end1, + Dictionary_it begin2, Dictionary_it end2, + Filtration_value filtration) { + if (begin1 == end1 || begin2 == end2) + return; // ----->> + while (true) { + if (begin1->first == begin2->first) { + intersection.push_back( + std::pair<Vertex_handle, Node>( + begin1->first, + Node(NULL, maximum(begin1->second.filtration(), begin2->second.filtration(), filtration)))); + ++begin1; + ++begin2; + if (begin1 == end1 || begin2 == end2) + return; // ----->> + } else { + if (begin1->first < begin2->first) { + ++begin1; + if (begin1 == end1) + return; + } else { + ++begin2; + if (begin2 == end2) + return; // ----->> + } + } + } + } + /** Maximum over 3 values.*/ + Filtration_value maximum(Filtration_value a, Filtration_value b, + Filtration_value c) { + Filtration_value max = (a < b) ? b : a; + return ((max < c) ? c : max); + } + + public: + /** \brief Write the hasse diagram of the simplicial complex in os. + * + * Each row in the file correspond to a simplex. A line is written: + * dim idx_1 ... idx_k fil where dim is the dimension of the simplex, + * idx_1 ... idx_k are the row index (starting from 0) of the simplices of the boundary + * of the simplex, and fil is its filtration value. */ + void print_hasse(std::ostream& os) { + os << num_simplices() << " " << std::endl; + for (auto sh : filtration_simplex_range(Indexing_tag())) { + os << dimension(sh) << " "; + for (auto b_sh : boundary_simplex_range(sh)) { + os << key(b_sh) << " "; + } + os << filtration(sh) << " \n"; + } + } + + private: Vertex_handle null_vertex_; -/** \brief Upper bound on the filtration values of the simplices.*/ - Filtration_value threshold_ ; -/** \brief Total number of simplices in the complex, without the empty simplex.*/ - size_t num_simplices_ ; -/** \brief Set of simplex tree Nodes representing the vertices.*/ - Siblings root_ ; -/** \brief Simplices ordered according to a filtration.*/ - std::vector< Simplex_handle > filtration_vect_ ; -/** \brief Upper bound on the dimension of the simplicial complex.*/ - int dimension_ ; + /** \brief Upper bound on the filtration values of the simplices.*/ + Filtration_value threshold_; + /** \brief Total number of simplices in the complex, without the empty simplex.*/ + unsigned int num_simplices_; + /** \brief Set of simplex tree Nodes representing the vertices.*/ + Siblings root_; + /** \brief Simplices ordered according to a filtration.*/ + std::vector<Simplex_handle> filtration_vect_; + /** \brief Upper bound on the dimension of the simplicial complex.*/ + int dimension_; }; // Print a Simplex_tree in os. -template< typename T1, typename T2, typename T3 > -std::ostream& operator<< ( std::ostream & os - , Simplex_tree< T1, T2, T3 > & st ) -{ - for(auto sh : st.filtration_simplex_range()) - { - os << st.dimension(sh) << " "; - for(auto v : st.simplex_vertex_range(sh)) - { os << v << " ";} - os << st.filtration(sh) << "\n";// TODO: VR - why adding the key ?? not read ?? << " " << st.key(sh) << " \n"; +template<typename T1, typename T2, typename T3> +std::ostream& operator<<(std::ostream & os, Simplex_tree<T1, T2, T3> & st) { + for (auto sh : st.filtration_simplex_range()) { + os << st.dimension(sh) << " "; + for (auto v : st.simplex_vertex_range(sh)) { + os << v << " "; } + os << st.filtration(sh) << "\n"; // TODO(VR): why adding the key ?? not read ?? << " " << st.key(sh) << " \n"; + } return os; } -template< typename T1, typename T2, typename T3 > -std::istream& operator>> ( std::istream & is - , Simplex_tree< T1, T2, T3 > & st ) -{ - //assert(st.num_simplices() == 0); - - std::vector< typename Simplex_tree<T1,T2,T3>::Vertex_handle > simplex; - typename Simplex_tree<T1,T2,T3>::Filtration_value fil; - typename Simplex_tree<T1,T2,T3>::Filtration_value max_fil = 0; - int max_dim = -1; - size_t num_simplices = 0; - while(read_simplex( is, simplex, fil )) //read all simplices in the file as a list of vertices - { +template<typename T1, typename T2, typename T3> +std::istream& operator>>(std::istream & is, Simplex_tree<T1, T2, T3> & st) { + // assert(st.num_simplices() == 0); + + std::vector<typename Simplex_tree<T1, T2, T3>::Vertex_handle> simplex; + typename Simplex_tree<T1, T2, T3>::Filtration_value fil; + typename Simplex_tree<T1, T2, T3>::Filtration_value max_fil = 0; + int max_dim = -1; + size_t num_simplices = 0; + while (read_simplex(is, simplex, fil)) { // read all simplices in the file as a list of vertices ++num_simplices; - int dim = (int)simplex.size() - 1; // Warning : simplex_size needs to be casted in int - Can be 0 - if(max_dim < dim) { max_dim = dim;} - if(max_fil < fil) { max_fil = fil;} - st.insert(simplex,fil); //insert every simplex in the simplex tree + int dim = static_cast<int>(simplex.size() - 1); // Warning : simplex_size needs to be casted in int - Can be 0 + if (max_dim < dim) { + max_dim = dim; + } + if (max_fil < fil) { + max_fil = fil; + } + st.insert(simplex, fil); // insert every simplex in the simplex tree simplex.clear(); } st.set_num_simplices(num_simplices); @@ -803,8 +842,7 @@ std::istream& operator>> ( std::istream & is return is; } -/** @} */ //end defgroup simplex_tree - -} // namespace GUDHI +/** @} */ // end defgroup simplex_tree +} // namespace Gudhi -#endif // GUDHI_FLAG_SIMPLEX_TREE_H +#endif // SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_ diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h index 6441a80a..06462c88 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h @@ -1,287 +1,306 @@ - /* 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): Clément Maria - * - * 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 SIMPLEX_TREE_ITERATORS_H -#define SIMPLEX_TREE_ITERATORS_H - -#include "boost/iterator/iterator_facade.hpp" - -namespace Gudhi{ +/* 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): Clément Maria + * + * 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 SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_ITERATORS_H_ +#define SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_ITERATORS_H_ + +#include <boost/iterator/iterator_facade.hpp> + +#include <vector> + +namespace Gudhi { /* \addtogroup simplex_tree - * Iterators and range types for the Simplex_tree. - * @{ - */ + * Iterators and range types for the Simplex_tree. + * @{ + */ /* \brief Iterator over the vertices of a simplex - * in a SimplexTree. - * - * Forward iterator, 'value_type' is SimplexTree::Vertex_handle.*/ -template < class SimplexTree > -class Simplex_tree_simplex_vertex_iterator -: public boost::iterator_facade < Simplex_tree_simplex_vertex_iterator < SimplexTree > - , typename SimplexTree::Vertex_handle const - , boost::forward_traversal_tag - , typename SimplexTree::Vertex_handle const - > -{ -public: + * in a SimplexTree. + * + * Forward iterator, 'value_type' is SimplexTree::Vertex_handle.*/ +template<class SimplexTree> +class Simplex_tree_simplex_vertex_iterator : public boost::iterator_facade< + Simplex_tree_simplex_vertex_iterator<SimplexTree>, + typename SimplexTree::Vertex_handle const, boost::forward_traversal_tag, + typename SimplexTree::Vertex_handle const> { + public: typedef typename SimplexTree::Simplex_handle Simplex_handle; - typedef typename SimplexTree::Siblings Siblings; - typedef typename SimplexTree::Vertex_handle Vertex_handle; + typedef typename SimplexTree::Siblings Siblings; + typedef typename SimplexTree::Vertex_handle Vertex_handle; - Simplex_tree_simplex_vertex_iterator (SimplexTree * st) : //any end() iterator - sib_(NULL), v_(st->null_vertex()) {} + explicit Simplex_tree_simplex_vertex_iterator(SimplexTree * st) + : // any end() iterator + sib_(NULL), + v_(st->null_vertex()) { + } - Simplex_tree_simplex_vertex_iterator( SimplexTree * st, - Simplex_handle sh) : - sib_(st->self_siblings(sh)), - v_(sh->first) {} + Simplex_tree_simplex_vertex_iterator(SimplexTree * st, Simplex_handle sh) + : sib_(st->self_siblings(sh)), + v_(sh->first) { + } -private: + private: friend class boost::iterator_core_access; - bool equal (Simplex_tree_simplex_vertex_iterator const &other) const - { return sib_ == other.sib_ && v_ == other.v_; } + bool equal(Simplex_tree_simplex_vertex_iterator const &other) const { + return sib_ == other.sib_ && v_ == other.v_; + } - Vertex_handle const& dereference() const { return v_; } + Vertex_handle const& dereference() const { + return v_; + } - void increment () { v_ = sib_->parent(); sib_ = sib_->oncles();} + void increment() { + v_ = sib_->parent(); + sib_ = sib_->oncles(); + } - Siblings * sib_; - Vertex_handle v_; + Siblings * sib_; + Vertex_handle v_; }; /*---------------------------------------------------------------------------*/ /* \brief Iterator over the simplices of the boundary of a - * simplex. - * - * Forward iterator, value_type is SimplexTree::Simplex_handle.*/ -template < class SimplexTree > -class Simplex_tree_boundary_simplex_iterator -: public boost::iterator_facade < Simplex_tree_boundary_simplex_iterator< SimplexTree > - , typename SimplexTree::Simplex_handle const - , boost::forward_traversal_tag - > -{ -public: - typedef typename SimplexTree::Simplex_handle Simplex_handle; - typedef typename SimplexTree::Vertex_handle Vertex_handle; - typedef typename SimplexTree::Siblings Siblings; + * simplex. + * + * Forward iterator, value_type is SimplexTree::Simplex_handle.*/ +template<class SimplexTree> +class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< + Simplex_tree_boundary_simplex_iterator<SimplexTree>, + typename SimplexTree::Simplex_handle const, boost::forward_traversal_tag> { + public: + typedef typename SimplexTree::Simplex_handle Simplex_handle; + typedef typename SimplexTree::Vertex_handle Vertex_handle; + typedef typename SimplexTree::Siblings Siblings; // any end() iterator - Simplex_tree_boundary_simplex_iterator(SimplexTree * st) : - last_(st->null_vertex()), sib_(NULL) {} - - Simplex_tree_boundary_simplex_iterator ( SimplexTree * st, - Simplex_handle sh ) : - suffix_(), st_(st) - { - last_ = sh->first; + explicit Simplex_tree_boundary_simplex_iterator(SimplexTree * st) + : last_(st->null_vertex()), + sib_(NULL) { + } + + Simplex_tree_boundary_simplex_iterator(SimplexTree * st, Simplex_handle sh) + : suffix_(), + st_(st) { + last_ = sh->first; Siblings * sib = st->self_siblings(sh); - next_ = sib->parent(); - sib_ = sib->oncles(); /* \todo check if NULL*/ - if(sib_ != NULL) { sh_ = sib_->find(next_); } - else { last_ = st->null_vertex(); } // vertex: == end() + next_ = sib->parent(); + sib_ = sib->oncles(); /* \todo check if NULL*/ + if (sib_ != NULL) { + sh_ = sib_->find(next_); + } else { + last_ = st->null_vertex(); + } // vertex: == end() } -private: + private: friend class boost::iterator_core_access; // valid when iterating along the SAME boundary. - bool equal (Simplex_tree_boundary_simplex_iterator const& other) const - { return (sib_ == other.sib_ && last_ == other.last_);} + bool equal(Simplex_tree_boundary_simplex_iterator const& other) const { + return (sib_ == other.sib_ && last_ == other.last_); + } + + Simplex_handle const& dereference() const { + return sh_; + } - Simplex_handle const& dereference () const { return sh_; } + void increment() { + if (sib_ == NULL) { + last_ = st_->null_vertex(); + return; + } - void increment() - { - if(sib_ == NULL) { last_ = st_->null_vertex(); return; } - Siblings * for_sib = sib_; - for(typename std::vector< Vertex_handle >::reverse_iterator rit = suffix_.rbegin(); - rit != suffix_.rend(); ++rit) - { + for (typename std::vector<Vertex_handle>::reverse_iterator rit = suffix_ + .rbegin(); rit != suffix_.rend(); ++rit) { sh_ = for_sib->find(*rit); - for_sib = sh_->second.children(); - } - sh_ = for_sib->find(last_); //sh_ points to the right simplex now + for_sib = sh_->second.children(); + } + sh_ = for_sib->find(last_); // sh_ points to the right simplex now suffix_.push_back(next_); next_ = sib_->parent(); sib_ = sib_->oncles(); } - Vertex_handle last_ ; //last vertex of the simplex - Vertex_handle next_ ; //next vertex to push in suffix_ - std::vector< Vertex_handle > suffix_ ; - Siblings * sib_ ; //where the next search will start from - Simplex_handle sh_ ; //current Simplex_handle in the boundary - SimplexTree * st_ ; //simplex containing the simplicial complex + Vertex_handle last_; // last vertex of the simplex + Vertex_handle next_; // next vertex to push in suffix_ + std::vector<Vertex_handle> suffix_; + Siblings * sib_; // where the next search will start from + Simplex_handle sh_; // current Simplex_handle in the boundary + SimplexTree * st_; // simplex containing the simplicial complex }; /*---------------------------------------------------------------------------*/ /* \brief Iterator over the simplices of a simplicial complex. - * - * Forward iterator, value_type is SimplexTree::Simplex_handle.*/ -template < class SimplexTree > -class Simplex_tree_complex_simplex_iterator -: public boost::iterator_facade < Simplex_tree_complex_simplex_iterator< SimplexTree > - , typename SimplexTree::Simplex_handle const - , boost::forward_traversal_tag - > -{ -public: + * + * Forward iterator, value_type is SimplexTree::Simplex_handle.*/ +template<class SimplexTree> +class Simplex_tree_complex_simplex_iterator : public boost::iterator_facade< + Simplex_tree_complex_simplex_iterator<SimplexTree>, + typename SimplexTree::Simplex_handle const, boost::forward_traversal_tag> { + public: typedef typename SimplexTree::Simplex_handle Simplex_handle; - typedef typename SimplexTree::Siblings Siblings; - typedef typename SimplexTree::Vertex_handle Vertex_handle; - -//any end() iterator - Simplex_tree_complex_simplex_iterator() : st_(NULL) {} - - Simplex_tree_complex_simplex_iterator(SimplexTree * st) : - st_(st) - { - if(st == NULL || st->root() == NULL || st->root()->members().empty()) { st_ = NULL; } - else - { - sh_ = st->root()->members().begin(); - sib_ = st->root(); - while(st->has_children(sh_)) - { sib_ = sh_->second.children(); - sh_ = sib_->members().begin(); } + typedef typename SimplexTree::Siblings Siblings; + typedef typename SimplexTree::Vertex_handle Vertex_handle; + +// any end() iterator + Simplex_tree_complex_simplex_iterator() + : st_(NULL) { + } + + explicit Simplex_tree_complex_simplex_iterator(SimplexTree * st) + : st_(st) { + if (st == NULL || st->root() == NULL || st->root()->members().empty()) { + st_ = NULL; + } else { + sh_ = st->root()->members().begin(); + sib_ = st->root(); + while (st->has_children(sh_)) { + sib_ = sh_->second.children(); + sh_ = sib_->members().begin(); + } } } -private: + private: friend class boost::iterator_core_access; // valid when iterating along the SAME boundary. - bool equal (Simplex_tree_complex_simplex_iterator const& other) const - { - if(other.st_ == NULL) { return (st_ == NULL); } - if(st_ == NULL) { return false; } + bool equal(Simplex_tree_complex_simplex_iterator const& other) const { + if (other.st_ == NULL) { + return (st_ == NULL); + } + if (st_ == NULL) { + return false; + } return (&(sh_->second) == &(other.sh_->second)); } - Simplex_handle const& dereference () const { return sh_; } + Simplex_handle const& dereference() const { + return sh_; + } // Depth first traversal. - void increment () - { + void increment() { ++sh_; - if(sh_ == sib_->members().end()) - { - if(sib_->oncles() == NULL) { st_ = NULL; return; } //reach the end + if (sh_ == sib_->members().end()) { + if (sib_->oncles() == NULL) { + st_ = NULL; + return; + } // reach the end sh_ = sib_->oncles()->members().find(sib_->parent()); - sib_ = sib_->oncles(); - return; + sib_ = sib_->oncles(); + return; } - while(st_->has_children(sh_)) - { - sib_ = sh_->second.children(); - sh_ = sib_->members().begin(); + while (st_->has_children(sh_)) { + sib_ = sh_->second.children(); + sh_ = sib_->members().begin(); } } - Simplex_handle sh_; - Siblings * sib_; - SimplexTree * st_; + Simplex_handle sh_; + Siblings * sib_; + SimplexTree * st_; }; /* \brief Iterator over the simplices of the skeleton of a given - * dimension of the simplicial complex. - * - * Forward iterator, value_type is SimplexTree::Simplex_handle.*/ -template < class SimplexTree > -class Simplex_tree_skeleton_simplex_iterator -: public boost::iterator_facade < Simplex_tree_skeleton_simplex_iterator< SimplexTree > - , typename SimplexTree::Simplex_handle const - , boost::forward_traversal_tag - > -{ -public: + * dimension of the simplicial complex. + * + * Forward iterator, value_type is SimplexTree::Simplex_handle.*/ +template<class SimplexTree> +class Simplex_tree_skeleton_simplex_iterator : public boost::iterator_facade< + Simplex_tree_skeleton_simplex_iterator<SimplexTree>, + typename SimplexTree::Simplex_handle const, boost::forward_traversal_tag> { + public: typedef typename SimplexTree::Simplex_handle Simplex_handle; - typedef typename SimplexTree::Siblings Siblings; - typedef typename SimplexTree::Vertex_handle Vertex_handle; - -//any end() iterator - Simplex_tree_skeleton_simplex_iterator() : st_(NULL) {} - - Simplex_tree_skeleton_simplex_iterator( SimplexTree * st - , int dim_skel ) - : st_(st) - , dim_skel_(dim_skel) - , curr_dim_(0) - { - if(st == NULL || st->root() == NULL || st->root()->members().empty()) { st_ = NULL; } - else - { - sh_ = st->root()->members().begin(); - sib_ = st->root(); - while(st->has_children(sh_) && curr_dim_ < dim_skel_) - { sib_ = sh_->second.children(); + typedef typename SimplexTree::Siblings Siblings; + typedef typename SimplexTree::Vertex_handle Vertex_handle; + +// any end() iterator + Simplex_tree_skeleton_simplex_iterator() + : st_(NULL) { + } + + Simplex_tree_skeleton_simplex_iterator(SimplexTree * st, int dim_skel) + : st_(st), + dim_skel_(dim_skel), + curr_dim_(0) { + if (st == NULL || st->root() == NULL || st->root()->members().empty()) { + st_ = NULL; + } else { + sh_ = st->root()->members().begin(); + sib_ = st->root(); + while (st->has_children(sh_) && curr_dim_ < dim_skel_) { + sib_ = sh_->second.children(); sh_ = sib_->members().begin(); - ++curr_dim_; } + ++curr_dim_; + } } } -private: + private: friend class boost::iterator_core_access; // valid when iterating along the SAME boundary. - bool equal (Simplex_tree_skeleton_simplex_iterator const& other) const - { - if(other.st_ == NULL) { return (st_ == NULL); } - if(st_ == NULL) { return false; } + bool equal(Simplex_tree_skeleton_simplex_iterator const& other) const { + if (other.st_ == NULL) { + return (st_ == NULL); + } + if (st_ == NULL) { + return false; + } return (&(sh_->second) == &(other.sh_->second)); } - Simplex_handle const& dereference () const { return sh_; } + Simplex_handle const& dereference() const { + return sh_; + } // Depth first traversal of the skeleton. - void increment () - { + void increment() { ++sh_; - if(sh_ == sib_->members().end()) - { - if(sib_->oncles() == NULL) { st_ = NULL; return; } //reach the end + if (sh_ == sib_->members().end()) { + if (sib_->oncles() == NULL) { + st_ = NULL; + return; + } // reach the end sh_ = sib_->oncles()->members().find(sib_->parent()); sib_ = sib_->oncles(); - --curr_dim_; - return; + --curr_dim_; + return; } - while(st_->has_children(sh_) && curr_dim_ < dim_skel_) - { - sib_ = sh_->second.children(); - sh_ = sib_->members().begin(); - ++curr_dim_; + while (st_->has_children(sh_) && curr_dim_ < dim_skel_) { + sib_ = sh_->second.children(); + sh_ = sib_->members().begin(); + ++curr_dim_; } } - Simplex_handle sh_; - Siblings * sib_; - SimplexTree * st_; - int dim_skel_; - int curr_dim_; + Simplex_handle sh_; + Siblings * sib_; + SimplexTree * st_; + int dim_skel_; + int curr_dim_; }; -/* @} */ //end addtogroup simplex_tree - -} // namespace GUDHI +/* @} */ // end addtogroup simplex_tree +} // namespace Gudhi -#endif // SIMPLEX_TREE_ITERATORS_H +#endif // SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_ITERATORS_H_ diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h index 91f4ef75..1f1a34cc 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h @@ -1,37 +1,36 @@ - /* 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): Clément Maria - * - * 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 GUDHI_SIMPLEX_TREE_NODE_EXPLICIT_STORAGE_H -#define GUDHI_SIMPLEX_TREE_NODE_EXPLICIT_STORAGE_H +/* 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): Clément Maria + * + * 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 SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_NODE_EXPLICIT_STORAGE_H_ +#define SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_NODE_EXPLICIT_STORAGE_H_ #include <vector> -#include <iostream> -namespace Gudhi{ +namespace Gudhi { /* \addtogroup simplex_tree - * Represents a node of a Simplex_tree. - * @{ - */ + * Represents a node of a Simplex_tree. + * @{ + */ /* * \brief Node of a simplex tree with filtration value @@ -39,66 +38,66 @@ namespace Gudhi{ * * It stores explicitely its own filtration value and its own Simplex_key. */ -template < class SimplexTree > +template<class SimplexTree> class Simplex_tree_node_explicit_storage { - public: -// friend SimplexTree; - - typedef typename SimplexTree::Siblings Siblings; + public: + typedef typename SimplexTree::Siblings Siblings; typedef typename SimplexTree::Filtration_value Filtration_value; - typedef typename SimplexTree::Simplex_key Simplex_key; + typedef typename SimplexTree::Simplex_key Simplex_key; - //private: - //friend class Simplex_tree; // Default constructor. - Simplex_tree_node_explicit_storage() : - children_(NULL), - simplex_key_(-1), - filtration_(0) {} + Simplex_tree_node_explicit_storage() + : children_(NULL), + simplex_key_(-1), + filtration_(0) { + } Simplex_tree_node_explicit_storage(Siblings * sib, - Filtration_value filtration) : - children_(sib), - simplex_key_(-1), - filtration_(filtration) {} + Filtration_value filtration) + : children_(sib), + simplex_key_(-1), + filtration_(filtration) { + } + void assign_key(Simplex_key key) { + simplex_key_ = key; + } - void assign_key(Simplex_key key) { simplex_key_ = key; } - - /* - * Return true if the node has children, - * false otherwise. - */ - //bool has_children(Vertex label) - //{ //if(children_ == NULL) return false; //for root simplices - // return (children_->parent() == label);} /* * Assign a children to the node */ - void assign_children (Siblings * children) { children_ = children; } + void assign_children(Siblings * children) { + children_ = children; + } /* * */ - void assign_filtration(double filtration_value) { filtration_ = filtration_value; } - - Filtration_value filtration() { return filtration_; } - - /* Careful -> has_children() must be true*/ - Siblings * children() { return children_; } - - Simplex_key key() { return simplex_key_; } - -private: - Siblings * children_; - + void assign_filtration(double filtration_value) { + filtration_ = filtration_value; + } + + Filtration_value filtration() { + return filtration_; + } + + /* Careful -> children_ can be NULL*/ + Siblings * children() { + return children_; + } + + Simplex_key key() { + return simplex_key_; + } + + private: + Siblings * children_; + // Data attached to simplex, explicit storage - Simplex_key simplex_key_; - Filtration_value filtration_; //value in the filtration - + Simplex_key simplex_key_; + Filtration_value filtration_; // value in the filtration }; -/* @} */ //end addtogroup simplex_tree - -} // namespace GUDHI +/* @} */ // end addtogroup simplex_tree +} // namespace Gudhi -#endif // GUDHI_SIMPLEX_TREE_NODE_EXPLICIT_STORAGE_H +#endif // SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_NODE_EXPLICIT_STORAGE_H_ diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h index 6db0d80a..977fafa1 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h @@ -1,87 +1,87 @@ - /* 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): Clément Maria - * - * 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 GUDHI_SIMPLEX_TREE_SIBLINGS -#define GUDHI_SIMPLEX_TREE_SIBLINGS +/* 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): Clément Maria + * + * 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 SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_SIBLINGS_H_ +#define SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_SIBLINGS_H_ #include "boost/container/flat_map.hpp" #include "Simplex_tree_node_explicit_storage.h" -namespace Gudhi{ +#include <utility> +#include <vector> + +namespace Gudhi { /* \addtogroup simplex_tree - * Represents a set of node of a Simplex_tree that share the same parent. - * @{ - */ + * Represents a set of node of a Simplex_tree that share the same parent. + * @{ + */ /* \brief Data structure to store a set of nodes in a SimplexTree sharing * the same parent node.*/ -template < class SimplexTree - , class MapContainer > +template<class SimplexTree, class MapContainer> class Simplex_tree_siblings { // private: // friend SimplexTree; -public: - template < class T > friend class Simplex_tree_simplex_vertex_iterator ; - template < class T > friend class Simplex_tree_boundary_simplex_iterator; - template < class T > friend class Simplex_tree_complex_simplex_iterator ; - template < class T > friend class Simplex_tree_skeleton_simplex_iterator; - - typedef typename SimplexTree::Vertex_handle Vertex_handle; - typedef typename SimplexTree::Filtration_value Filtration_value; - typedef typename SimplexTree::Node Node; - typedef MapContainer Dictionary; - typedef typename MapContainer::iterator Dictionary_it; + public: + template<class T> friend class Simplex_tree_simplex_vertex_iterator; + template<class T> friend class Simplex_tree_boundary_simplex_iterator; + template<class T> friend class Simplex_tree_complex_simplex_iterator; + template<class T> friend class Simplex_tree_skeleton_simplex_iterator; + + typedef typename SimplexTree::Vertex_handle Vertex_handle; + typedef typename SimplexTree::Filtration_value Filtration_value; + typedef typename SimplexTree::Node Node; + typedef MapContainer Dictionary; + typedef typename MapContainer::iterator Dictionary_it; /* Default constructor.*/ - Simplex_tree_siblings() - : oncles_(NULL) - , parent_(-1) - , members_() {} - + Simplex_tree_siblings() + : oncles_(NULL), + parent_(-1), + members_() { + } + /* Constructor with values.*/ - Simplex_tree_siblings(Simplex_tree_siblings * oncles, - Vertex_handle parent ) - : oncles_(oncles) - , parent_(parent) - , members_() {} - + Simplex_tree_siblings(Simplex_tree_siblings * oncles, Vertex_handle parent) + : oncles_(oncles), + parent_(parent), + members_() { + } + /* \brief Constructor with initialized set of members. - * - * 'members' must be sorted and unique.*/ - Simplex_tree_siblings(Simplex_tree_siblings * oncles, - Vertex_handle parent, - std::vector< std::pair< Vertex_handle, Node > > & members) - : oncles_ ( oncles ) - , parent_ ( parent ) - , members_ ( boost::container::ordered_unique_range - , members.begin() - , members.end() ) - { - for(auto map_it = members_.begin(); - map_it != members_.end(); map_it++) - { map_it->second.assign_children(this); } + * + * 'members' must be sorted and unique.*/ + Simplex_tree_siblings(Simplex_tree_siblings * oncles, Vertex_handle parent, + const std::vector<std::pair<Vertex_handle, Node> > & members) + : oncles_(oncles), + parent_(parent), + members_(boost::container::ordered_unique_range, members.begin(), + members.end()) { + for (auto map_it = members_.begin(); map_it != members_.end(); map_it++) { + map_it->second.assign_children(this); + } } - + /* * \brief Inserts a Node in the set of siblings nodes. * @@ -89,41 +89,45 @@ public: * between input filtration_value and the value already * present in the node. */ - void insert ( Vertex_handle v - , Filtration_value filtration_value ) - { + void insert(Vertex_handle v, Filtration_value filtration_value) { typename Dictionary::iterator sh = members_.find(v); - if(sh != members_.end() && sh->second.filtration() > filtration_value) - { sh->second.assign_filtration(filtration_value); - return; } - if(sh == members_.end()) - { members_.insert(std::pair< Vertex_handle, Node >( v, Node(this,filtration_value) )); - return; } + if (sh != members_.end() && sh->second.filtration() > filtration_value) { + sh->second.assign_filtration(filtration_value); + return; + } + if (sh == members_.end()) { + members_.insert( + std::pair<Vertex_handle, Node>(v, Node(this, filtration_value))); + return; + } } - typename Dictionary::iterator find( Vertex_handle v ) - { return members_.find(v); } - - Simplex_tree_siblings * oncles() - { return oncles_; } - - Vertex_handle parent() - { return parent_; } - - Dictionary & members() - { return members_; } - - size_t size() { return members_.size(); } - - - Simplex_tree_siblings * oncles_ ; - Vertex_handle parent_ ; - Dictionary members_ ; - -}; + typename Dictionary::iterator find(Vertex_handle v) { + return members_.find(v); + } + + Simplex_tree_siblings * oncles() { + return oncles_; + } + + Vertex_handle parent() { + return parent_; + } -/* @} */ //end addtogroup simplex_tree + Dictionary & members() { + return members_; + } + + size_t size() { + return members_.size(); + } + + Simplex_tree_siblings * oncles_; + Vertex_handle parent_; + Dictionary members_; +}; -} // namespace GUDHI +/* @} */ // end addtogroup simplex_tree +} // namespace Gudhi -#endif // GUDHI_SIMPLEX_TREE_SIBLINGS +#endif // SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_SIMPLEX_TREE_SIBLINGS_H_ diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/indexing_tag.h b/src/Simplex_tree/include/gudhi/Simplex_tree/indexing_tag.h index 165d166d..69ffa44b 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/indexing_tag.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/indexing_tag.h @@ -1,34 +1,39 @@ - /* 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): Clément Maria - * - * 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/>. - */ +/* 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): Clément Maria + * + * 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/>. + */ -namespace Gudhi{ +#ifndef SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_INDEXING_TAG_H_ +#define SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_INDEXING_TAG_H_ + +namespace Gudhi { /** \brief Tag for a linear ordering of simplices. - * - * \implements IndexingTag - */ - struct linear_indexing_tag {}; + * + * \implements IndexingTag + */ +struct linear_indexing_tag { +}; /* \brief Tag for a zigzag ordering of simplices. */ // struct zigzag_indexing_tag {}; +} // namespace Gudhi -} // namespace GUDHI +#endif // SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_INDEXING_TAG_H_ diff --git a/src/Simplex_tree/test/CMakeLists.txt b/src/Simplex_tree/test/CMakeLists.txt index a15ac04e..02ef9d8b 100644 --- a/src/Simplex_tree/test/CMakeLists.txt +++ b/src/Simplex_tree/test/CMakeLists.txt @@ -1,11 +1,21 @@ cmake_minimum_required(VERSION 2.6) project(GUDHITestSimplexTree) -# NEED TO FIND BOOST NEEDED COMPONENTS TO LINK WITH -find_package(Boost 1.45.0 COMPONENTS system unit_test_framework) -message("Boost_lib = ${Boost_LIBRARIES}") +if(NOT MSVC) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +endif() -include_directories(${Boost_INCLUDE_DIRS}) -add_executable ( TestSimplexTree UnitTestSimplexTree.cpp ) -target_link_libraries(TestSimplexTree ${Boost_LIBRARIES}) +add_executable ( simplex_tree_unit_test simplex_tree_unit_test.cpp ) +target_link_libraries(simplex_tree_unit_test ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +# Unitary tests +add_test(simplex_tree_unit_test ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_unit_test) + +if (LCOV_PATH) + # Lcov code coverage of unitary test + add_test(src/Simplex_tree/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Simplex_tree) +endif() + +cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Simplex_tree/include/gudhi") diff --git a/src/Simplex_tree/test/README b/src/Simplex_tree/test/README index 3d6981ff..620bcd5f 100644 --- a/src/Simplex_tree/test/README +++ b/src/Simplex_tree/test/README @@ -7,6 +7,6 @@ make To launch with details: *********************** -./TestSimplexTree --report_level=detailed --log_level=all +./simplex_tree_unit_test --report_level=detailed --log_level=all ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Simplex_tree/test/UnitTestSimplexTree.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index a3671f56..b44f95a8 100644 --- a/src/Simplex_tree/test/UnitTestSimplexTree.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -166,7 +166,7 @@ void set_and_test_simplex_tree_dim_fil(typeST& simplexTree, int vectorSize, cons simplexTree.set_filtration(max_fil); std::cout << " set_and_test_simplex_tree_dim_fil - max_fil=" << max_fil << std::endl; } - int nb_simplices = simplexTree.num_simplices() + 1; + unsigned int nb_simplices = simplexTree.num_simplices() + 1; simplexTree.set_num_simplices(nb_simplices); BOOST_CHECK( simplexTree.dimension() == dim_max ); @@ -347,11 +347,12 @@ BOOST_AUTO_TEST_CASE( simplex_tree_insertion ) BOOST_CHECK( st.dimension() == dim_max ); BOOST_CHECK( AreAlmostTheSame(st.filtration(), max_fil) ); - // 1 - // o - // /X\ - // o---o---o - // 2 0 3 + /* Inserted simplex: */ + /* 1 */ + /* o */ + /* /X\ */ + /* o---o---o */ + /* 2 0 3 */ // [0.1] 0 // [0.1] 1 diff --git a/src/Skeleton_blocker/example/CMakeLists.txt b/src/Skeleton_blocker/example/CMakeLists.txt index 0aed1e7a..c8e976cd 100644 --- a/src/Skeleton_blocker/example/CMakeLists.txt +++ b/src/Skeleton_blocker/example/CMakeLists.txt @@ -1,7 +1,7 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIskbl) -add_executable ( SkeletonBlockerIteration Skeleton_blocker_iteration.cpp ) - -target_link_libraries( SkeletonBlockerIteration ${Boost_TIMER_LIBRARY}) +add_executable(SkeletonBlockerIteration Skeleton_blocker_iteration.cpp) +target_link_libraries(SkeletonBlockerIteration ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) +add_test(SkeletonBlockerIteration ${CMAKE_CURRENT_BINARY_DIR}/SkeletonBlockerIteration) diff --git a/src/Skeleton_blocker/test/CMakeLists.txt b/src/Skeleton_blocker/test/CMakeLists.txt index d66dcf64..c341dcee 100644 --- a/src/Skeleton_blocker/test/CMakeLists.txt +++ b/src/Skeleton_blocker/test/CMakeLists.txt @@ -1,7 +1,21 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIskbl) -add_executable ( TestSkeletonBlockerComplex TestSkeletonBlockerComplex.cpp ) -add_executable ( TestSimplifiable TestSimplifiable.cpp ) +if(NOT MSVC) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +endif() +add_executable(TestSkeletonBlockerComplex TestSkeletonBlockerComplex.cpp) +add_executable(TestSimplifiable TestSimplifiable.cpp) +add_test(TestSkeletonBlockerComplex ${CMAKE_CURRENT_BINARY_DIR}/TestSkeletonBlockerComplex) +add_test(TestSimplifiable ${CMAKE_CURRENT_BINARY_DIR}/TestSimplifiable) + +if (LCOV_PATH) + # Lcov code coverage of unitary test + add_test(src/Skeleton_blocker/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Skeleton_blocker) +endif() + +cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Skeleton_blocker/include/gudhi") |