Merge of cpplint branch

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@357 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 8b6af51684401f4c45c953142aae90d6f45e6778

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")