Merged with trunk and removed unnecessary files

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

Former-commit-id: d0ec52d222d22c102e9fe57590882cd0024c82d5

14 files changed, 884 insertions, 777 deletions

diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt
index 9487cce6..50d10025 100644
--- a/src/Persistent_cohomology/example/CMakeLists.txt
+++ b/src/Persistent_cohomology/example/CMakeLists.txt
@@ -2,39 +2,47 @@ cmake_minimum_required(VERSION 2.6)
 project(GUDHIExPersCohom)
 
 # problem with Visual Studio link on Boost program_options
-if (NOT MSVC)
-   add_executable(persistence_from_simple_simplex_tree persistence_from_simple_simplex_tree.cpp)
-   target_link_libraries(persistence_from_simple_simplex_tree ${Boost_SYSTEM_LIBRARY})
+add_definitions( -DBOOST_ALL_NO_LIB )
+add_definitions( -DBOOST_ALL_DYN_LINK )
 
-   add_executable(rips_persistence rips_persistence.cpp)
-   target_link_libraries(rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
-   
-   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(persistence_from_file persistence_from_file.cpp)
-   target_link_libraries(persistence_from_file ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
-   add_test(persistence_from_file_3_2_0 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 2 -m 0)
-   add_test(persistence_from_file_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 3 -m 100)
-   
-   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_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES})
-      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(plain_homology plain_homology.cpp)
+target_link_libraries(plain_homology ${Boost_SYSTEM_LIBRARY})
+add_test(plain_homology ${CMAKE_CURRENT_BINARY_DIR}/plain_homology)
+
+add_executable(persistence_from_simple_simplex_tree persistence_from_simple_simplex_tree.cpp)
+target_link_libraries(persistence_from_simple_simplex_tree ${Boost_SYSTEM_LIBRARY})
+add_test(persistence_from_simple_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_simple_simplex_tree 1 0)
+
+add_executable(rips_persistence rips_persistence.cpp)
+target_link_libraries(rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
    
-      add_executable ( performance_rips_persistence performance_rips_persistence.cpp )
-      target_link_libraries(performance_rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES})
+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)
    
-      if(CGAL_FOUND)
-         # uncomment to display debug traces
-         # add_definitions(-DDEBUG_TRACES)
-         add_executable(alpha_shapes_persistence alpha_shapes_persistence.cpp)
-         target_link_libraries(alpha_shapes_persistence ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY})
-         add_test(alpha_shapes_persistence_2_0_5 ${CMAKE_CURRENT_BINARY_DIR}/alpha_shapes_persistence ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 2 0.5)
-         #add_test(alpha_shapes_persistence_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/alpha_shapes_persistence ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 3 -m 100)
-      endif()
+add_executable(persistence_from_file persistence_from_file.cpp)
+target_link_libraries(persistence_from_file ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+add_test(persistence_from_file_3_2_0 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 2 -m 0)
+add_test(persistence_from_file_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 3 -m 100)
    
+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_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES})
+   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_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES})
+
+   if(CGAL_FOUND)
+     if (CMAKE_BUILD_TYPE MATCHES Debug)
+       # For programs to be more verbose
+       add_definitions(-DDEBUG_TRACES)
+     endif()
+     add_executable(alpha_shapes_persistence alpha_shapes_persistence.cpp)
+     target_link_libraries(alpha_shapes_persistence ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY})
+     add_test(alpha_shapes_persistence_2_0_5 ${CMAKE_CURRENT_BINARY_DIR}/alpha_shapes_persistence ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 2 0.5)
+     #add_test(alpha_shapes_persistence_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/alpha_shapes_persistence ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 3 -m 100)
    endif()
+   
 endif()
diff --git a/src/Persistent_cohomology/example/alpha_shapes_persistence.cpp b/src/Persistent_cohomology/example/alpha_shapes_persistence.cpp
index 3a50c84c..6d5eebcf 100644
--- a/src/Persistent_cohomology/example/alpha_shapes_persistence.cpp
+++ b/src/Persistent_cohomology/example/alpha_shapes_persistence.cpp
@@ -1,24 +1,29 @@
- /*    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):       Vincent Rouvreau
-  *
-  *    Copyright (C) 2014  INRIA Saclay (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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014  INRIA Saclay (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 <gudhi/graph_simplicial_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <boost/variant.hpp>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Delaunay_triangulation_3.h>
@@ -27,36 +32,37 @@
 
 #include <fstream>
 #include <cmath>
-
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.h"
-#include <boost/variant.hpp>
+#include <string>
+#include <tuple>
+#include <map>
+#include <utility>
+#include <list>
+#include <vector>
 
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
 
 // Alpha_shape_3 templates type definitions
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef CGAL::Alpha_shape_vertex_base_3<Kernel>             Vb;
-typedef CGAL::Alpha_shape_cell_base_3<Kernel>               Fb;
-typedef CGAL::Triangulation_data_structure_3<Vb,Fb>         Tds;
-typedef CGAL::Delaunay_triangulation_3<Kernel,Tds>          Triangulation_3;
-typedef CGAL::Alpha_shape_3<Triangulation_3>                Alpha_shape_3;
+typedef CGAL::Alpha_shape_vertex_base_3<Kernel> Vb;
+typedef CGAL::Alpha_shape_cell_base_3<Kernel> Fb;
+typedef CGAL::Triangulation_data_structure_3<Vb, Fb> Tds;
+typedef CGAL::Delaunay_triangulation_3<Kernel, Tds> Triangulation_3;
+typedef CGAL::Alpha_shape_3<Triangulation_3> Alpha_shape_3;
 
 // From file type definition
-typedef Kernel::Point_3                                     Point_3;
+typedef Kernel::Point_3 Point_3;
 
 // filtration with alpha values needed type definition
 typedef Alpha_shape_3::FT Alpha_value_type;
-typedef CGAL::Object      Object;
+typedef CGAL::Object Object;
 typedef CGAL::Dispatch_output_iterator<
-  CGAL::cpp11::tuple<Object, Alpha_value_type>,
-  CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, std::back_insert_iterator< std::vector<Alpha_value_type> >
-                     > > Dispatch;
-typedef Alpha_shape_3::Cell_handle   Cell_handle;
-typedef Alpha_shape_3::Facet         Facet;
-typedef Alpha_shape_3::Edge          Edge_3;
+CGAL::cpp11::tuple<Object, Alpha_value_type>,
+CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >,
+    std::back_insert_iterator< std::vector<Alpha_value_type> > > > Dispatch;
+typedef Alpha_shape_3::Cell_handle Cell_handle;
+typedef Alpha_shape_3::Facet Facet;
+typedef Alpha_shape_3::Edge Edge_3;
 typedef std::list<Alpha_shape_3::Vertex_handle> Vertex_list;
 
 // gudhi type definition
@@ -65,70 +71,60 @@ typedef std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex > Alpha_shape
 typedef std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex> Alpha_shape_simplex_tree_pair;
 typedef std::vector< Simplex_tree_vertex > Simplex_tree_vector_vertex;
 
-//#define DEBUG_TRACES
-
-Vertex_list from (const Cell_handle& ch)
-{
+Vertex_list from(const Cell_handle& ch) {
   Vertex_list the_list;
-  for (auto i = 0; i < 4; i++)
-  {
+  for (auto i = 0; i < 4; i++) {
 #ifdef DEBUG_TRACES
     std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
     the_list.push_back(ch->vertex(i));
   }
   return the_list;
 }
-Vertex_list from (const Facet& fct)
-{
+
+Vertex_list from(const Facet& fct) {
   Vertex_list the_list;
-  for (auto i = 0; i < 4; i++)
-  {
-    if (fct.second != i)
-    {
+  for (auto i = 0; i < 4; i++) {
+    if (fct.second != i) {
 #ifdef DEBUG_TRACES
       std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
       the_list.push_back(fct.first->vertex(i));
     }
   }
   return the_list;
 }
-Vertex_list from (const Edge_3& edg)
-{
+
+Vertex_list from(const Edge_3& edg) {
   Vertex_list the_list;
-  for (auto i = 0; i < 4; i++)
-  {
-    if ((edg.second == i) ||(edg.third == i))
-    {
+  for (auto i = 0; i < 4; i++) {
+    if ((edg.second == i) || (edg.third == i)) {
 #ifdef DEBUG_TRACES
       std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
       the_list.push_back(edg.first->vertex(i));
     }
   }
   return the_list;
 }
-Vertex_list from (const Alpha_shape_3::Vertex_handle& vh)
-{
+
+Vertex_list from(const Alpha_shape_3::Vertex_handle& vh) {
   Vertex_list the_list;
 #ifdef DEBUG_TRACES
   std::cout << "from vertex=" << vh->point() << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
   the_list.push_back(vh);
   return the_list;
 }
 
-void usage(char * const progName)
-{
-  std::cerr << "Usage: " << progName << " path_to_file_graph coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n";
-  exit(-1); // ----- >>
+void usage(char * const progName) {
+  std::cerr << "Usage: " << progName <<
+      " path_to_file_graph coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n";
+  exit(-1);
 }
 
-int main (int argc, char * const argv[])
-{
-
-  int coeff_field_characteristic=0;
+int main(int argc, char * const argv[]) {
+  int coeff_field_characteristic = 0;
   int returnedScanValue = sscanf(argv[2], "%d", &coeff_field_characteristic);
   if ((returnedScanValue == EOF) || (coeff_field_characteristic <= 0)) {
     std::cerr << "Error: " << argv[2] << " is not correct\n";
@@ -149,102 +145,96 @@ int main (int argc, char * const argv[])
   }
 
   // Read points from file
-  std::string filegraph   = argv[1];
+  std::string filegraph = argv[1];
   std::list<Point_3> lp;
   std::ifstream is(filegraph.c_str());
   int n;
   is >> n;
 #ifdef DEBUG_TRACES
   std::cout << "Reading " << n << " points " << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
   Point_3 p;
-  for( ; n>0 ; n--)    {
+  for (; n > 0; n--) {
     is >> p;
     lp.push_back(p);
   }
 
   // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
-  Alpha_shape_3 as(lp.begin(),lp.end(),0,Alpha_shape_3::GENERAL);
+  Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL);
 #ifdef DEBUG_TRACES
   std::cout << "Alpha shape computed in GENERAL mode" << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
 
   // filtration with alpha values from alpha shape
   std::vector<Object> the_objects;
   std::vector<Alpha_value_type> the_alpha_values;
 
-  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>( std::back_inserter(the_objects), std::back_inserter(the_alpha_values));
+  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
+                                                                  std::back_inserter(the_alpha_values));
 
   as.filtration_with_alpha_values(disp);
 #ifdef DEBUG_TRACES
   std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
 
   Alpha_shape_3::size_type count_vertices = 0;
-  Alpha_shape_3::size_type count_edges    = 0;
-  Alpha_shape_3::size_type count_facets   = 0;
-  Alpha_shape_3::size_type count_cells    = 0;
+  Alpha_shape_3::size_type count_edges = 0;
+  Alpha_shape_3::size_type count_facets = 0;
+  Alpha_shape_3::size_type count_cells = 0;
 
   // Loop on objects vector
   Vertex_list vertex_list;
   Simplex_tree<> simplex_tree;
   Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
   std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
-  int dim_max=0;
-  Filtration_value filtration_max=0.0;
-  for(auto object_iterator: the_objects)
-  {
+  int dim_max = 0;
+  Filtration_value filtration_max = 0.0;
+  for (auto object_iterator : the_objects) {
     // Retrieve Alpha shape vertex list from object
-    if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator))
-    {
+    if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
       vertex_list = from(*cell);
       count_cells++;
       if (dim_max < 3) {
-        dim_max=3; // Cell is of dim 3
+        // Cell is of dim 3
+        dim_max = 3;
       }
-    }
-    else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator))
-    {
+    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
       vertex_list = from(*facet);
       count_facets++;
       if (dim_max < 2) {
-        dim_max=2; // Facet is of dim 2
+        // Facet is of dim 2
+        dim_max = 2;
       }
-    }
-    else if (const Edge_3* edge = CGAL::object_cast<Edge_3>(&object_iterator))
-    {
+    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
       vertex_list = from(*edge);
       count_edges++;
       if (dim_max < 1) {
-        dim_max=1; // Edge_3 is of dim 1
+        // Edge_3 is of dim 1
+        dim_max = 1;
       }
-    }
-    else if (const Alpha_shape_3::Vertex_handle* vertex = CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator))
-    {
+    } else if (const Alpha_shape_3::Vertex_handle * vertex =
+               CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) {
       count_vertices++;
       vertex_list = from(*vertex);
     }
     // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
     Simplex_tree_vector_vertex the_simplex_tree;
-    for (auto the_alpha_shape_vertex:vertex_list)
-    {
+    for (auto the_alpha_shape_vertex : vertex_list) {
       Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
-      if (the_map_iterator == map_cgal_simplex_tree.end())
-      {
+      if (the_map_iterator == map_cgal_simplex_tree.end()) {
         // alpha shape not found
         Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
 #ifdef DEBUG_TRACES
         std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
         the_simplex_tree.push_back(vertex);
-        map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex,vertex));
-      } else
-      {
+        map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex));
+      } else {
         // alpha shape found
         Simplex_tree_vertex vertex = the_map_iterator->second;
 #ifdef DEBUG_TRACES
         std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
         the_simplex_tree.push_back(vertex);
       }
     }
@@ -252,7 +242,7 @@ int main (int argc, char * const argv[])
     Filtration_value filtr = std::sqrt(*the_alpha_value_iterator);
 #ifdef DEBUG_TRACES
     std::cout << "filtration = " << filtr << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
     if (filtr > filtration_max) {
       filtration_max = filtr;
     }
@@ -263,14 +253,13 @@ int main (int argc, char * const argv[])
       std::cout << "This shall not happen" << std::endl;
   }
   simplex_tree.set_filtration(filtration_max);
-  simplex_tree.set_num_simplices(count_vertices + count_edges + count_facets + count_cells);
   simplex_tree.set_dimension(dim_max);
 
 #ifdef DEBUG_TRACES
   std::cout << "vertices \t\t" << count_vertices << std::endl;
-  std::cout << "edges \t\t"    << count_edges << std::endl;
-  std::cout << "facets \t\t"   << count_facets << std::endl;
-  std::cout << "cells \t\t"    << count_cells << std::endl;
+  std::cout << "edges \t\t" << count_edges << std::endl;
+  std::cout << "facets \t\t" << count_facets << std::endl;
+  std::cout << "cells \t\t" << count_cells << std::endl;
 
 
   std::cout << "Information of the Simplex Tree: " << std::endl;
@@ -278,23 +267,25 @@ int main (int argc, char * const argv[])
   std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
   std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
   std::cout << "  filtration = " << simplex_tree.filtration() << std::endl << std::endl;
-#endif // DEBUG_TRACES
+#endif  // DEBUG_TRACES
 
 #ifdef DEBUG_TRACES
   std::cout << "Iterator on vertices: " << std::endl;
-  for( auto vertex : simplex_tree.complex_vertex_range() )
-  { std::cout << vertex << " "; }
-#endif // DEBUG_TRACES
+  for (auto vertex : simplex_tree.complex_vertex_range()) {
+    std::cout << vertex << " ";
+  }
+#endif  // DEBUG_TRACES
 
   // Sort the simplices in the order of the filtration
   simplex_tree.initialize_filtration();
 
   std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
   // Compute the persistence diagram of the complex
-  Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh( simplex_tree );
-  pcoh.init_coefficients( coeff_field_characteristic ); //initializes the coefficient field for homology
+  Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
 
-  pcoh.compute_persistent_cohomology( min_persistence );
+  pcoh.compute_persistent_cohomology(min_persistence);
 
   pcoh.output_diagram();
 
diff --git a/src/Persistent_cohomology/example/performance_rips_persistence.cpp b/src/Persistent_cohomology/example/performance_rips_persistence.cpp
index 077c2b07..fc48d6b1 100644
--- a/src/Persistent_cohomology/example/performance_rips_persistence.cpp
+++ b/src/Persistent_cohomology/example/performance_rips_persistence.cpp
@@ -1,49 +1,51 @@
- /*    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/>.
-  */
-
-#include "gudhi/reader_utils.h"
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.h"
-#include "gudhi/Persistent_cohomology/Multi_field.h"
-#include "gudhi/Hasse_complex.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/>.
+ */
+
+#include <gudhi/reader_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Persistent_cohomology/Multi_field.h>
+#include <gudhi/Hasse_complex.h>
 
 #include <chrono>
+#include <string>
+#include <vector>
 
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
 
 /* Compute the persistent homology of the complex cpx with coefficients in Z/pZ. */
 template< typename FilteredComplex>
-void timing_persistence( FilteredComplex & cpx
-                       , int               p );
+void timing_persistence(FilteredComplex & cpx
+                        , int p);
 
 /* Compute multi-field persistent homology of the complex cpx with coefficients in
  * Z/rZ for all prime number r in [p;q].*/
 template< typename FilteredComplex>
-void timing_persistence( FilteredComplex & cpx
-                       , int               p 
-                       , int               q );
+void timing_persistence(FilteredComplex & cpx
+                        , int p
+                        , int q);
 
 /* Timings for the computation of persistent homology with different 
  * representations of a Rips complex and different coefficient fields. The 
@@ -59,111 +61,154 @@ void timing_persistence( FilteredComplex & cpx
  * We present also timings for the computation of multi-field persistent 
  * homology in all fields Z/rZ for r prime between 2 and 1223.
  */
-int main (int argc, char * argv[]) 
-{
+int main(int argc, char * argv[]) {
   std::chrono::time_point<std::chrono::system_clock> start, end;
-  int enlapsed_sec;
+  int elapsed_sec;
+  {
 
-  std::string      filepoints      = "../examples/Kl.txt";
-  Filtration_value threshold       = 0.3;
-  int              dim_max         = 3;
-  int              p               = 2;
-  int              q               = 1223;
+  std::string filepoints = "../../../data/points/Kl.txt";
+  Filtration_value threshold = 0.27;
+  int dim_max = 3;
+  int p = 2;
+  int q = 1223;
 
-// Extract the points from the file filepoints
+  // Extract the points from the file filepoints
   typedef std::vector<double> Point_t;
   std::vector< Point_t > points;
-  read_points( filepoints, points );
+  read_points(filepoints, points);
 
-// Compute the proximity graph of the points
+  // Compute the proximity graph of the points
   start = std::chrono::system_clock::now();
-  Graph_t prox_graph = compute_proximity_graph( points, threshold
-                                              , euclidean_distance<Point_t> );
+  Graph_t prox_graph = compute_proximity_graph(points, threshold
+                                               , euclidean_distance<Point_t>);
   end = std::chrono::system_clock::now();
-  enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count();
-  std::cout << "Compute Rips graph in " << enlapsed_sec << " sec.\n";
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n";
 
-// Construct the Rips complex in a Simplex Tree
-  Simplex_tree<> st;        
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree<> st;
   start = std::chrono::system_clock::now();
 
-  st.insert_graph(prox_graph); // insert the proximity graph in the simplex tree
-  st.expansion( dim_max ); // expand the graph until dimension dim_max
+  // insert the proximity graph in the simplex tree
+  st.insert_graph(prox_graph);
+  // expand the graph until dimension dim_max
+  st.expansion(dim_max);
 
   end = std::chrono::system_clock::now();
-  enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count();
-  std::cout << "Compute Rips complex in " << enlapsed_sec << " sec.\n";
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Compute Rips complex in " << elapsed_sec << " ms.\n";
   std::cout << "  - dimension           = " << st.dimension() << std::endl;
   std::cout << "  - number of simplices = " << st.num_simplices() << std::endl;
 
-// Sort the simplices in the order of the filtration
+  // Sort the simplices in the order of the filtration
   start = std::chrono::system_clock::now();
   st.initialize_filtration();
   end = std::chrono::system_clock::now();
-  enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count();
-  std::cout << "Order the simplices of the filtration in " << enlapsed_sec << " sec.\n";
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Order the simplices of the filtration in " << elapsed_sec << " ms.\n";
+
+  // Copy the keys inside the simplices
+  start = std::chrono::system_clock::now();
+  {
+    int count = 0;
+    for (auto sh : st.filtration_simplex_range())
+      st.assign_key(sh, count++);
+  }
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Copied the keys inside the simplices in " << elapsed_sec << " ms.\n";
 
-// Convert the simplex tree into a hasse diagram
+  // Convert the simplex tree into a hasse diagram
   start = std::chrono::system_clock::now();
   Hasse_complex<> hcpx(st);
   end = std::chrono::system_clock::now();
-  enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count();
-  std::cout << "Convert the simplex tree into a Hasse diagram in " << enlapsed_sec << " sec.\n";
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Convert the simplex tree into a Hasse diagram in " << elapsed_sec << " ms.\n";
 
 
   std::cout << "Timings when using a simplex tree: \n";
-  timing_persistence(st,p);
-  timing_persistence(st,q);
-  timing_persistence(st,p,q);
+  timing_persistence(st, p);
+  timing_persistence(st, q);
+  timing_persistence(st, p, q);
 
   std::cout << "Timings when using a Hasse complex: \n";
-  timing_persistence(hcpx,p);
-  timing_persistence(hcpx,q);
-  timing_persistence(hcpx,p,q);
+  timing_persistence(hcpx, p);
+  timing_persistence(hcpx, q);
+  timing_persistence(hcpx, p, q);
 
+  start = std::chrono::system_clock::now();
+  }
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "Running the complex destructors in " << elapsed_sec << " ms.\n";
   return 0;
 }
 
-
 template< typename FilteredComplex>
 void
-timing_persistence( FilteredComplex & cpx
-                  , int               p )
-{
+timing_persistence(FilteredComplex & cpx
+                   , int p) {
   std::chrono::time_point<std::chrono::system_clock> start, end;
-  int enlapsed_sec;
+  int elapsed_sec;
+  {
+  start = std::chrono::system_clock::now();
+  Persistent_cohomology< FilteredComplex, Field_Zp > pcoh(cpx);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Initialize pcoh in " << elapsed_sec << " ms.\n";
+  // initializes the coefficient field for homology
+  start = std::chrono::system_clock::now();
+  pcoh.init_coefficients(p);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Initialize the coefficient field in " << elapsed_sec << " ms.\n";
 
-  Persistent_cohomology< FilteredComplex, Field_Zp > pcoh (cpx);
-  pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-  
   start = std::chrono::system_clock::now();
- 
-  pcoh.compute_persistent_cohomology( INFINITY );
- 
+
+  pcoh.compute_persistent_cohomology(INFINITY);
+
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Compute persistent homology in Z/" << p << "Z in " << elapsed_sec << " ms.\n";
+  start = std::chrono::system_clock::now();
+  }
   end = std::chrono::system_clock::now();
-  enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count();
-  std::cout << "  Compute persistent homology in Z/"<<p<<"Z in " << enlapsed_sec << " sec.\n";
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Run the persistence destructors in " << elapsed_sec << " ms.\n";
 }
 
 template< typename FilteredComplex>
 void
-timing_persistence( FilteredComplex & cpx
-                  , int               p 
-                  , int               q )
-{
+timing_persistence(FilteredComplex & cpx
+                   , int p
+                   , int q) {
   std::chrono::time_point<std::chrono::system_clock> start, end;
-  int enlapsed_sec;
-
-  Persistent_cohomology< FilteredComplex, Multi_field > pcoh (cpx);
-  pcoh.init_coefficients( p, q ); //initilizes the coefficient field for homology
+  int elapsed_sec;
+  {
+  start = std::chrono::system_clock::now();
+  Persistent_cohomology< FilteredComplex, Multi_field > pcoh(cpx);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Initialize pcoh in " << elapsed_sec << " ms.\n";
+  // initializes the coefficient field for homology
+  start = std::chrono::system_clock::now();
+  pcoh.init_coefficients(p, q);
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Initialize the coefficient field in " << elapsed_sec << " ms.\n";
   // compute persistent homology, disgarding persistent features of life shorter than min_persistence
 
   start = std::chrono::system_clock::now();
 
-  pcoh.compute_persistent_cohomology( INFINITY );
+  pcoh.compute_persistent_cohomology(INFINITY);
 
   end = std::chrono::system_clock::now();
-  enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
   std::cout << "  Compute multi-field persistent homology in all coefficient fields Z/pZ "
-            << "with p in ["<<p<<";"<<q<<"] in " << enlapsed_sec << " sec.\n";
+      << "with p in [" << p << ";" << q << "] in " << elapsed_sec << " ms.\n";
+  start = std::chrono::system_clock::now();
+  }
+  end = std::chrono::system_clock::now();
+  elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+  std::cout << "  Run the persistence destructors in " << elapsed_sec << " ms.\n";
 }
diff --git a/src/Persistent_cohomology/example/persistence_from_file.cpp b/src/Persistent_cohomology/example/persistence_from_file.cpp
index e886aea7..8eb8d0f3 100644
--- a/src/Persistent_cohomology/example/persistence_from_file.cpp
+++ b/src/Persistent_cohomology/example/persistence_from_file.cpp
@@ -1,141 +1,144 @@
- /*    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):       Vincent Rouvreau
-  *
-  *    Copyright (C) 2014  INRIA Saclay (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 "gudhi/reader_utils.h"
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014  INRIA Saclay (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 <gudhi/reader_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
 
 #include <boost/program_options.hpp>
 
+#include <string>
+
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
 
-typedef int        Vertex_handle;
-typedef double     Filtration_value;
-
-void program_options( int argc, char *   argv[] 
-                    , std::string      & simplex_tree_file
-                    , std::string      & output_file
-                    , int              & p
-                    , Filtration_value & min_persistence );
-
-int main (int argc, char * argv[]) 
-{
-	std::string      simplex_tree_file;
-	std::string      output_file  ;
-	int              p         ;
-	Filtration_value min_persistence;
-
-	program_options(argc,argv,simplex_tree_file,output_file,p,min_persistence);
-
-	std::cout << "Simplex_tree from file=" << simplex_tree_file.c_str() << " - output_file=" << output_file.c_str() << std::endl;
-	std::cout << "     - p=" << p << " - min_persistence=" << min_persistence << std::endl;
-
-	// Construct the Rips complex in a Simplex Tree
-	Simplex_tree<> simplex_tree;
-
-	std::ifstream simplex_tree_stream(simplex_tree_file);
-	simplex_tree_stream >> simplex_tree;
-
-	std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices" << std::endl;
-	std::cout << "   - dimension " << simplex_tree.dimension() << "   - filtration " << simplex_tree.filtration() << std::endl;
-
-	/*
-	std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl;
-	for( auto f_simplex : simplex_tree.filtration_simplex_range() )
-	{ std::cout << "   " << "[" << simplex_tree.filtration(f_simplex) << "] ";
-	for( auto vertex : simplex_tree.simplex_vertex_range(f_simplex) )
-	{ std::cout << vertex << " "; }
-	std::cout << std::endl;
-	}*/
-
-	// Sort the simplices in the order of the filtration
-	simplex_tree.initialize_filtration();
-
-	// Compute the persistence diagram of the complex
-	Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh( simplex_tree );
-	pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-
-	pcoh.compute_persistent_cohomology( min_persistence );
-
-	// Output the diagram in output_file
-	if(output_file.empty()) {  pcoh.output_diagram(); }
-	else {
-		std::ofstream out(output_file);
-		pcoh.output_diagram(out);
-		out.close(); }
+typedef int Vertex_handle;
+typedef double Filtration_value;
+
+void program_options(int argc, char * argv[]
+                     , std::string & simplex_tree_file
+                     , std::string & output_file
+                     , int & p
+                     , Filtration_value & min_persistence);
+
+int main(int argc, char * argv[]) {
+  std::string simplex_tree_file;
+  std::string output_file;
+  int p;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, simplex_tree_file, output_file, p, min_persistence);
+
+  std::cout << "Simplex_tree from file=" << simplex_tree_file.c_str() << " - output_file=" << output_file.c_str()
+      << std::endl;
+  std::cout << "     - p=" << p << " - min_persistence=" << min_persistence << std::endl;
+
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree<> simplex_tree;
+
+  std::ifstream simplex_tree_stream(simplex_tree_file);
+  simplex_tree_stream >> simplex_tree;
+
+  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices" << std::endl;
+  std::cout << "   - dimension " << simplex_tree.dimension() << "   - filtration " << simplex_tree.filtration()
+      << std::endl;
+
+  /*
+  std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl;
+  for( auto f_simplex : simplex_tree.filtration_simplex_range() )
+  { std::cout << "   " << "[" << simplex_tree.filtration(f_simplex) << "] ";
+  for( auto vertex : simplex_tree.simplex_vertex_range(f_simplex) )
+  { std::cout << vertex << " "; }
+  std::cout << std::endl;
+  }*/
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(p);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in output_file
+  if (output_file.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(output_file);
+    pcoh.output_diagram(out);
+    out.close();
+  }
 
-	return 0;
+  return 0;
 }
 
-
-
-void program_options( int argc, char *   argv[] 
-                    , std::string      & simplex_tree_file
-                    , std::string      & output_file
-                    , int              & p
-                    , Filtration_value & min_persistence )
-{
+void program_options(int argc, char * argv[]
+                     , std::string & simplex_tree_file
+                     , std::string & output_file
+                     , int & p
+                     , Filtration_value & min_persistence) {
   namespace po = boost::program_options;
-  po::options_description     hidden("Hidden options");
+  po::options_description hidden("Hidden options");
   hidden.add_options()
-        ("input-file",  po::value<std::string>(&simplex_tree_file),
-          "Name of file containing a simplex set. Format is one simplex per line (cf. reader_utils.h - read_simplex):   Dim1 X11 X12 ... X1d Fil1  ");
-    
+      ("input-file", po::value<std::string>(&simplex_tree_file),
+       "Name of file containing a simplex set. Format is one simplex per line (cf. reader_utils.h - read_simplex): Dim1 X11 X12 ... X1d Fil1  ");
+
   po::options_description visible("Allowed options", 100);
   visible.add_options()
-        ("help,h",          "produce help message")
-        ("output-file,o", po::value<std::string>(&output_file)->default_value(std::string()),
-                   "Name of file in which the persistence diagram is written. Default print in std::cout")
-        ("field-charac,p", po::value<int>(&p)->default_value(11),                  
-                  "Characteristic p of the coefficient field Z/pZ for computing homology.")
-        ("min-persistence,m", po::value<Filtration_value>(&min_persistence), 
-                  "Minimal lifetime of homology feature to be recorded. Default is 0");
+      ("help,h", "produce help message")
+      ("output-file,o", po::value<std::string>(&output_file)->default_value(std::string()),
+       "Name of file in which the persistence diagram is written. Default print in std::cout")
+      ("field-charac,p", po::value<int>(&p)->default_value(11),
+       "Characteristic p of the coefficient field Z/pZ for computing homology.")
+      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
+       "Minimal lifetime of homology feature to be recorded. Default is 0");
 
   po::positional_options_description pos;
   pos.add("input-file", 1);
-  
-  po::options_description all; all.add(visible).add(hidden);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
 
   po::variables_map vm;
   po::store(po::command_line_parser(argc, argv).
-                options(all).positional(pos).run(), vm);
+            options(all).positional(pos).run(), vm);
   po::notify(vm);
 
-  if (vm.count("help") || !vm.count("input-file"))
-  { 
-      std::cout << std::endl;
-      std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-      std::cout << "of a Rips complex defined on a set of input points.\n \n";
-      std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-      std::cout << "   p   dim b d \n";
-      std::cout << "where dim is the dimension of the homological feature,\n";
-      std::cout << "b and d are respectively the birth and death of the feature and \n";
-      std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
-
-      std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-      std::cout << visible << std::endl; 
-      std::abort();
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a Rips complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
   }
 }
diff --git a/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp
index ba82e4e6..ba772f04 100644
--- a/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp
+++ b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp
@@ -20,11 +20,14 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+
 #include <iostream>
 #include <ctime>
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.h"
+#include <utility>
+#include <vector>
 
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
@@ -34,15 +37,19 @@ typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex;
 typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool;
 typedef Simplex_tree<> typeST;
 
-void usage(char * const progName)
-{
+void usage(char * const progName) {
   std::cerr << "Usage: " << progName << " coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n";
-  exit(-1); // ----- >>
+  exit(-1);
 }
 
-int main (int argc, char * const argv[])
-{
-  int coeff_field_characteristic=0;
+int main(int argc, char * const argv[]) {
+  // program args management
+  if (argc != 3) {
+    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
+    usage(argv[0]);
+  }
+
+  int coeff_field_characteristic = 0;
   int returnedScanValue = sscanf(argv[1], "%d", &coeff_field_characteristic);
   if ((returnedScanValue == EOF) || (coeff_field_characteristic <= 0)) {
     std::cerr << "Error: " << argv[1] << " is not correct\n";
@@ -56,101 +63,65 @@ int main (int argc, char * const argv[])
     usage(argv[0]);
   }
 
-  // program args management
-  if (argc != 3) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[0]);
-  }
-
   // TEST OF INSERTION
   std::cout << "********************************************************************" << std::endl;
   std::cout << "TEST OF INSERTION" << std::endl;
   typeST st;
 
   // ++ FIRST
-  std::cout << "   - INSERT (2,1,0)" << std::endl;
-  typeVectorVertex SimplexVector1;
-  SimplexVector1.push_back(2);
-  SimplexVector1.push_back(1);
-  SimplexVector1.push_back(0);
-  st.insert_simplex_and_subfaces ( SimplexVector1, 0.3);
+  std::cout << "   - INSERT (0,1,2)" << std::endl;
+  typeVectorVertex SimplexVector = {0, 1, 2};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.3);
 
   // ++ SECOND
   std::cout << "   - INSERT 3" << std::endl;
-  typeVectorVertex SimplexVector2;
-  SimplexVector2.push_back(3);
-  st.insert_simplex_and_subfaces ( SimplexVector2, 0.1);
+  SimplexVector = {3};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.1);
 
   // ++ THIRD
   std::cout << "   - INSERT (0,3)" << std::endl;
-  typeVectorVertex SimplexVector3;
-  SimplexVector3.push_back(3);
-  SimplexVector3.push_back(0);
-  st.insert_simplex_and_subfaces ( SimplexVector3, 0.2);
+  SimplexVector = {0, 3};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.2);
 
   // ++ FOURTH
-  std::cout << "   - INSERT (1,0) (already inserted)" << std::endl;
-  typeVectorVertex SimplexVector4;
-  SimplexVector4.push_back(1);
-  SimplexVector4.push_back(0);
-  st.insert_simplex_and_subfaces ( SimplexVector4, 0.2);
+  std::cout << "   - INSERT (0,1) (already inserted)" << std::endl;
+  SimplexVector = {0, 1};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.2);
 
   // ++ FIFTH
   std::cout << "   - INSERT (3,4,5)" << std::endl;
-  typeVectorVertex SimplexVector5;
-  SimplexVector5.push_back(3);
-  SimplexVector5.push_back(4);
-  SimplexVector5.push_back(5);
-  st.insert_simplex_and_subfaces ( SimplexVector5, 0.3);
+  SimplexVector = {3, 4, 5};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.3);
 
   // ++ SIXTH
   std::cout << "   - INSERT (0,1,6,7)" << std::endl;
-  typeVectorVertex SimplexVector6;
-  SimplexVector6.push_back(0);
-  SimplexVector6.push_back(1);
-  SimplexVector6.push_back(6);
-  SimplexVector6.push_back(7);
-  st.insert_simplex_and_subfaces ( SimplexVector6, 0.4);
+  SimplexVector = {0, 1, 6, 7};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.4);
 
   // ++ SEVENTH
   std::cout << "   - INSERT (4,5,8,9)" << std::endl;
-  typeVectorVertex SimplexVector7;
-  SimplexVector7.push_back(4);
-  SimplexVector7.push_back(5);
-  SimplexVector7.push_back(8);
-  SimplexVector7.push_back(9);
-  st.insert_simplex_and_subfaces ( SimplexVector7, 0.4);
+  SimplexVector = {4, 5, 8, 9};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.4);
 
   // ++ EIGHTH
   std::cout << "   - INSERT (9,10,11)" << std::endl;
-  typeVectorVertex SimplexVector8;
-  SimplexVector8.push_back(9);
-  SimplexVector8.push_back(10);
-  SimplexVector8.push_back(11);
-  st.insert_simplex_and_subfaces ( SimplexVector8, 0.3);
-  
+  SimplexVector = {9, 10, 11};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.3);
+
   // ++ NINETH
   std::cout << "   - INSERT (2,10,12)" << std::endl;
-  typeVectorVertex SimplexVector9;
-  SimplexVector9.push_back(2);
-  SimplexVector9.push_back(10);
-  SimplexVector9.push_back(12);
-  st.insert_simplex_and_subfaces ( SimplexVector9, 0.3);
-  
+  SimplexVector = {2, 10, 12};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.3);
+
   // ++ TENTH
   std::cout << "   - INSERT (11,6)" << std::endl;
-  typeVectorVertex SimplexVector10;
-  SimplexVector10.push_back(11);
-  SimplexVector10.push_back(6);
-  st.insert_simplex_and_subfaces ( SimplexVector10, 0.2);
+  SimplexVector = {6, 11};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.2);
 
   // ++ ELEVENTH
   std::cout << "   - INSERT (13,14,15)" << std::endl;
-  typeVectorVertex SimplexVector11;
-  SimplexVector11.push_back(13);
-  SimplexVector11.push_back(14);
-  SimplexVector11.push_back(15);
-  st.insert_simplex_and_subfaces ( SimplexVector11, 0.25);
+  SimplexVector = {13, 14, 15};
+  st.insert_simplex_and_subfaces(SimplexVector, 0.25);
 
   /* Inserted simplex:         */
   /*    1   6                  */
@@ -175,33 +146,31 @@ int main (int argc, char * const argv[])
   st.set_dimension(2);
   st.set_filtration(0.4);
 
-  std::cout << "The complex contains " << st.num_simplices() << " simplices - " << st.num_vertices() << " vertices " << std::endl;
+  std::cout << "The complex contains " << st.num_simplices() << " simplices - " << st.num_vertices() << " vertices "
+      << std::endl;
   std::cout << "   - dimension " << st.dimension() << "   - filtration " << st.filtration() << std::endl;
-  std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl;
+  std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:"
+      << std::endl;
   std::cout << "**************************************************************" << std::endl;
   std::cout << "strict graph G { " << std::endl;
 
-  for( auto f_simplex : st.filtration_simplex_range() )
-  {
+  for (auto f_simplex : st.filtration_simplex_range()) {
     std::cout << "   " << "[" << st.filtration(f_simplex) << "] ";
-    for( auto vertex : st.simplex_vertex_range(f_simplex) )
-    {
-      std::cout << (int)vertex << " -- ";
+    for (auto vertex : st.simplex_vertex_range(f_simplex)) {
+      std::cout << static_cast<int>(vertex) << " -- ";
     }
     std::cout << ";" << std::endl;
   }
 
   std::cout << "}" << std::endl;
-  //std::cout << "**************************************************************" << std::endl;
-  //st.print_hasse(std::cout);
   std::cout << "**************************************************************" << std::endl;
 
-
   // Compute the persistence diagram of the complex
   persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(st);
-  pcoh.init_coefficients( coeff_field_characteristic ); //initiliazes the coefficient field for homology
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
 
-  pcoh.compute_persistent_cohomology( min_persistence );
+  pcoh.compute_persistent_cohomology(min_persistence);
 
   // Output the diagram in filediag
   pcoh.output_diagram();
diff --git a/src/Persistent_cohomology/example/plain_homology.cpp b/src/Persistent_cohomology/example/plain_homology.cpp
new file mode 100644
index 00000000..0a692717
--- /dev/null
+++ b/src/Persistent_cohomology/example/plain_homology.cpp
@@ -0,0 +1,84 @@
+/*    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):       Marc Glisse
+ *
+ *    Copyright (C) 2015  INRIA Saclay - Ile-de-France (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 <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+
+#include <iostream>
+
+using namespace Gudhi;
+
+/* We could perfectly well use the default Simplex_tree<> (which uses
+ * Simplex_tree_options_full_featured), the following simply demonstrates
+ * how to save on storage by not storing a filtration value.  */
+
+struct MyOptions : Simplex_tree_options_full_featured {
+  // Implicitly use 0 as filtration value for all simplices
+  static const bool store_filtration = false;
+  // The persistence algorithm needs this
+  static const bool store_key = true;
+  // I have few vertices
+  typedef short Vertex_handle;
+};
+typedef Simplex_tree<MyOptions> ST;
+
+int main() {
+  ST st;
+
+  /* Complex to build. */
+  /*    1   3          */
+  /*    o---o          */
+  /*   /X\ /           */
+  /*  o---o   o        */
+  /*  2   0   4        */
+
+  const short triangle012[] = {0, 1, 2};
+  const short edge03[] = {0, 3};
+  const short edge13[] = {1, 3};
+  const short vertex4[] = {4};
+  st.insert_simplex_and_subfaces(triangle012);
+  st.insert_simplex_and_subfaces(edge03);
+  st.insert_simplex(edge13);
+  st.insert_simplex(vertex4);
+  // FIXME: Remove this line
+  st.set_dimension(2);
+
+  // Sort the simplices in the order of the filtration
+  st.initialize_filtration();
+
+  // Class for homology computation
+  persistent_cohomology::Persistent_cohomology<ST, persistent_cohomology::Field_Zp> pcoh(st);
+
+  // Initialize the coefficient field Z/2Z for homology
+  pcoh.init_coefficients(2);
+
+  // Compute the persistence diagram of the complex
+  pcoh.compute_persistent_cohomology();
+
+  // Print the result. The format is, on each line: 2 dim 0 inf
+  // where 2 represents the field, dim the dimension of the feature.
+  // 2  0 0 inf 
+  // 2  0 0 inf 
+  // 2  1 0 inf 
+  // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=1.
+  pcoh.output_diagram();
+}
diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp
index 297a8f98..5277bf7a 100644
--- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp
+++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp
@@ -1,153 +1,157 @@
- /*    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/>.
-  */
-
-#include "gudhi/reader_utils.h"
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.h"
-#include "gudhi/Persistent_cohomology/Multi_field.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/>.
+ */
+
+#include <gudhi/reader_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Persistent_cohomology/Multi_field.h>
 
 #include <boost/program_options.hpp>
 
+#include <string>
+#include <vector>
+
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
 
-typedef int    Vertex_handle;
+typedef int Vertex_handle;
 typedef double Filtration_value;
 
-void program_options( int argc, char *   argv[] 
-                    , std::string      & filepoints
-                    , std::string      & filediag
-                    , Filtration_value & threshold
-                    , int              & dim_max
-                    , int              & min_p
-                    , int              & max_p
-                    , Filtration_value & min_persistence );
-
-int main (int argc, char * argv[]) 
-{
-  std::string      filepoints;
-  std::string      filediag;
+void program_options(int argc, char * argv[]
+                     , std::string & filepoints
+                     , std::string & filediag
+                     , Filtration_value & threshold
+                     , int & dim_max
+                     , int & min_p
+                     , int & max_p
+                     , Filtration_value & min_persistence);
+
+int main(int argc, char * argv[]) {
+  std::string filepoints;
+  std::string filediag;
   Filtration_value threshold;
-  int              dim_max;
-  int              min_p;
-  int              max_p;
+  int dim_max;
+  int min_p;
+  int max_p;
   Filtration_value min_persistence;
 
-  program_options(argc,argv,filepoints,filediag,threshold,dim_max,min_p,max_p,min_persistence);
+  program_options(argc, argv, filepoints, filediag, threshold, dim_max, min_p, max_p, min_persistence);
 
-// Extract the points from the file filepoints
+  // Extract the points from the file filepoints
   typedef std::vector<double> Point_t;
   std::vector< Point_t > points;
-  read_points( filepoints, points );
+  read_points(filepoints, points);
 
-// Compute the proximity graph of the points
-  Graph_t prox_graph = compute_proximity_graph( points, threshold
-                                              , euclidean_distance<Point_t> );
+  // Compute the proximity graph of the points
+  Graph_t prox_graph = compute_proximity_graph(points, threshold
+                                               , euclidean_distance<Point_t>);
 
-// Construct the Rips complex in a Simplex Tree
-  Simplex_tree<> st;        
-  st.insert_graph(prox_graph); // insert the proximity graph in the simplex tree
-  st.expansion( dim_max ); // expand the graph until dimension dim_max
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree<> st;
+  // insert the proximity graph in the simplex tree
+  st.insert_graph(prox_graph);
+  // expand the graph until dimension dim_max
+  st.expansion(dim_max);
 
-// Sort the simplices in the order of the filtration
+  // 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_p, max_p ); //initilizes the coefficient field for homology
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology< Simplex_tree<>, Multi_field > pcoh(st);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(min_p, max_p);
   // compute persistent homology, disgarding persistent features of life shorter than min_persistence
-  pcoh.compute_persistent_cohomology( min_persistence );
-
-// Output the diagram in filediag
-  if(filediag.empty()) {  pcoh.output_diagram(); }
-  else {
-   std::ofstream out(filediag);
-   pcoh.output_diagram(out);
-   out.close(); }
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
 
   return 0;
 }
 
-
-
-void program_options( int argc, char *   argv[] 
-                    , std::string      & filepoints
-                    , std::string      & filediag
-                    , Filtration_value & threshold
-                    , int              & dim_max
-                    , int              & min_p
-                    , int              & max_p
-                    , Filtration_value & min_persistence )
-{
+void program_options(int argc, char * argv[]
+                     , std::string & filepoints
+                     , std::string & filediag
+                     , Filtration_value & threshold
+                     , int & dim_max
+                     , int & min_p
+                     , int & max_p
+                     , Filtration_value & min_persistence) {
   namespace po = boost::program_options;
-  po::options_description     hidden("Hidden options");
+  po::options_description hidden("Hidden options");
   hidden.add_options()
-        ("input-file",  po::value<std::string>(&filepoints), 
-          "Name of file containing a point set. Format is one point per line:   X1 ... Xd \n");
-    
+      ("input-file", po::value<std::string>(&filepoints),
+       "Name of file containing a point set. Format is one point per line:   X1 ... Xd \n");
+
   po::options_description visible("Allowed options");
   visible.add_options()
-        ("help,h",          "produce help message")
-        ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), 
-                   "Name of file in which the persistence diagram is written. Default print in std::cout")
-        ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0), 
-                   "Maximal length of an edge for the Rips complex construction.")
-        ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),        
-                  "Maximal dimension of the Rips complex we want to compute.")
-        ("min-field-charac,p", po::value<int>(&min_p)->default_value(2),                  
-                  "Minimal characteristic p of the coefficient field Z/pZ.")
-        ("max-field-charac,q", po::value<int>(&max_p)->default_value(1223),                  
-                  "Minimial characteristic q of the coefficient field Z/pZ.")
-        ("min-persistence,m", po::value<Filtration_value>(&min_persistence), 
-                  "Minimal lifetime of homology feature to be recorded. Default is 0");
+      ("help,h", "produce help message")
+      ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
+       "Name of file in which the persistence diagram is written. Default print in std::cout")
+      ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0),
+       "Maximal length of an edge for the Rips complex construction.")
+      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+       "Maximal dimension of the Rips complex we want to compute.")
+      ("min-field-charac,p", po::value<int>(&min_p)->default_value(2),
+       "Minimal characteristic p of the coefficient field Z/pZ.")
+      ("max-field-charac,q", po::value<int>(&max_p)->default_value(1223),
+       "Minimial characteristic q of the coefficient field Z/pZ.")
+      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
+       "Minimal lifetime of homology feature to be recorded. Default is 0");
 
   po::positional_options_description pos;
   pos.add("input-file", 1);
-  
-  po::options_description all; all.add(visible).add(hidden);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
 
   po::variables_map vm;
   po::store(po::command_line_parser(argc, argv).
-                options(all).positional(pos).run(), vm);
+            options(all).positional(pos).run(), vm);
   po::notify(vm);
 
-  if (vm.count("help") || !vm.count("input-file"))
-  { 
-      std::cout << std::endl;
-      std::cout << "Compute the persistent homology with various coefficient fields \n";
-      std::cout << "of a Rips complex defined on a set of input points. The coefficient \n";
-      std::cout << "fields are all the Z/rZ for a prime number r contained in the \n";
-      std::cout << "specified range [p,q]\n \n";
-      std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-      std::cout << "   p1*...*pr   dim b d \n";
-      std::cout << "where dim is the dimension of the homological feature,\n";
-      std::cout << "b and d are respectively the birth and death of the feature and \n";
-      std::cout << "p1*...*pr is the product of prime numbers pi such that the homology \n";
-      std::cout << "feature exists in homology with Z/piZ coefficients."<< std::endl << std::endl;
-
-      std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-      std::cout << visible << std::endl; 
-      std::abort();
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with various coefficient fields \n";
+    std::cout << "of a Rips complex defined on a set of input points. The coefficient \n";
+    std::cout << "fields are all the Z/rZ for a prime number r contained in the \n";
+    std::cout << "specified range [p,q]\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p1*...*pr   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p1*...*pr is the product of prime numbers pi such that the homology \n";
+    std::cout << "feature exists in homology with Z/piZ coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
   }
 }
diff --git a/src/Persistent_cohomology/example/rips_persistence.cpp b/src/Persistent_cohomology/example/rips_persistence.cpp
index 4253def9..9b1ef42f 100644
--- a/src/Persistent_cohomology/example/rips_persistence.cpp
+++ b/src/Persistent_cohomology/example/rips_persistence.cpp
@@ -1,147 +1,151 @@
- /*    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/>.
-  */
-
-#include "gudhi/reader_utils.h"
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include "gudhi/Persistent_cohomology.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/>.
+ */
+
+#include <gudhi/reader_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
 
 #include <boost/program_options.hpp>
 
+#include <string>
+#include <vector>
+
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
 
-typedef int        Vertex_handle;
-typedef double     Filtration_value;
-
-void program_options( int argc, char *   argv[] 
-                    , std::string      & filepoints
-                    , std::string      & filediag
-                    , Filtration_value & threshold
-                    , int              & dim_max
-                    , int              & p
-                    , Filtration_value & min_persistence );
-
-int main (int argc, char * argv[]) 
-{
-  std::string      filepoints;
-  std::string      filediag  ;
-  Filtration_value threshold ;
-  int              dim_max   ;
-  int              p         ;
+typedef int Vertex_handle;
+typedef double Filtration_value;
+
+void program_options(int argc, char * argv[]
+                     , std::string & filepoints
+                     , std::string & filediag
+                     , Filtration_value & threshold
+                     , int & dim_max
+                     , int & p
+                     , Filtration_value & min_persistence);
+
+int main(int argc, char * argv[]) {
+  std::string filepoints;
+  std::string filediag;
+  Filtration_value threshold;
+  int dim_max;
+  int p;
   Filtration_value min_persistence;
 
-  program_options(argc,argv,filepoints,filediag,threshold,dim_max,p,min_persistence);
+  program_options(argc, argv, filepoints, filediag, threshold, dim_max, p, min_persistence);
 
-// Extract the points from the file filepoints
+  // Extract the points from the file filepoints
   typedef std::vector<double> Point_t;
   std::vector< Point_t > points;
-  read_points( filepoints, points );
+  read_points(filepoints, points);
 
-// Compute the proximity graph of the points
-  Graph_t prox_graph = compute_proximity_graph( points, threshold
-                                              , euclidean_distance<Point_t> );
+  // Compute the proximity graph of the points
+  Graph_t prox_graph = compute_proximity_graph(points, threshold
+                                               , euclidean_distance<Point_t>);
 
-// Construct the Rips complex in a Simplex Tree
-  Simplex_tree<> st;        
-  st.insert_graph(prox_graph); // insert the proximity graph in the simplex tree
-  st.expansion( dim_max ); // expand the graph until dimension dim_max
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree<> st;
+  // insert the proximity graph in the simplex tree
+  st.insert_graph(prox_graph);
+  // expand the graph until dimension dim_max
+  st.expansion(dim_max);
 
   std::cout << "The complex contains " << st.num_simplices() << " simplices \n";
   std::cout << "   and has dimension " << st.dimension() << " \n";
 
-// Sort the simplices in the order of the filtration
+  // Sort the simplices in the order of the filtration
   st.initialize_filtration();
 
-// Compute the persistence diagram of the complex
+  // Compute the persistence diagram of the complex
   persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(st);
-  pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-  
-  pcoh.compute_persistent_cohomology( min_persistence );
-
-// Output the diagram in filediag
-  if(filediag.empty()) {  pcoh.output_diagram(); }
-  else {
-   std::ofstream out(filediag);
-   pcoh.output_diagram(out);
-   out.close(); }
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(p);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
 
   return 0;
 }
 
-
-
-void program_options( int argc, char *   argv[] 
-                    , std::string      & filepoints
-                    , std::string      & filediag
-                    , Filtration_value & threshold
-                    , int              & dim_max
-                    , int              & p
-                    , Filtration_value & min_persistence )
-{
+void program_options(int argc, char * argv[]
+                     , std::string & filepoints
+                     , std::string & filediag
+                     , Filtration_value & threshold
+                     , int & dim_max
+                     , int & p
+                     , Filtration_value & min_persistence) {
   namespace po = boost::program_options;
-  po::options_description     hidden("Hidden options");
+  po::options_description hidden("Hidden options");
   hidden.add_options()
-        ("input-file",  po::value<std::string>(&filepoints), 
-          "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
-    
+      ("input-file", po::value<std::string>(&filepoints),
+       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
+
   po::options_description visible("Allowed options", 100);
   visible.add_options()
-        ("help,h",          "produce help message")
-        ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), 
-                   "Name of file in which the persistence diagram is written. Default print in std::cout")
-        ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0), 
-                   "Maximal length of an edge for the Rips complex construction.")
-        ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),        
-                  "Maximal dimension of the Rips complex we want to compute.")
-        ("field-charac,p", po::value<int>(&p)->default_value(11),                  
-                  "Characteristic p of the coefficient field Z/pZ for computing homology.")
-        ("min-persistence,m", po::value<Filtration_value>(&min_persistence), 
-                  "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
+      ("help,h", "produce help message")
+      ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
+       "Name of file in which the persistence diagram is written. Default print in std::cout")
+      ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0),
+       "Maximal length of an edge for the Rips complex construction.")
+      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+       "Maximal dimension of the Rips complex we want to compute.")
+      ("field-charac,p", po::value<int>(&p)->default_value(11),
+       "Characteristic p of the coefficient field Z/pZ for computing homology.")
+      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
+       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
 
   po::positional_options_description pos;
   pos.add("input-file", 1);
-  
-  po::options_description all; all.add(visible).add(hidden);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
 
   po::variables_map vm;
   po::store(po::command_line_parser(argc, argv).
-                options(all).positional(pos).run(), vm);
+            options(all).positional(pos).run(), vm);
   po::notify(vm);
 
-  if (vm.count("help") || !vm.count("input-file"))
-  { 
-      std::cout << std::endl;
-      std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-      std::cout << "of a Rips complex defined on a set of input points.\n \n";
-      std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-      std::cout << "   p   dim b d \n";
-      std::cout << "where dim is the dimension of the homological feature,\n";
-      std::cout << "b and d are respectively the birth and death of the feature and \n";
-      std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
-
-      std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-      std::cout << visible << std::endl; 
-      std::abort();
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a Rips complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
   }
 }
diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
index b0d68f09..2a405830 100644
--- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
+++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
@@ -20,23 +20,28 @@
  *    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 PERSISTENT_COHOMOLOGY_H_
+#define PERSISTENT_COHOMOLOGY_H_
 
 #include <gudhi/Persistent_cohomology/Persistent_cohomology_column.h>
 #include <gudhi/Persistent_cohomology/Field_Zp.h>
+#include <gudhi/Simple_object_pool.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 <map>
 #include <utility>
 #include <list>
 #include <vector>
 #include <set>
+#include <fstream>  // std::ofstream
+#include <limits>  // for numeric_limits<>
+#include <tuple>
+#include <algorithm>
+#include <string>
 
 namespace Gudhi {
 
@@ -229,9 +234,10 @@ class Persistent_cohomology {
       : 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
+        num_simplices_(cpx_->num_simplices()),           // num_simplices save to avoid to call thrice the function
+        ds_rank_(num_simplices_),                        // union-find
+        ds_parent_(num_simplices_),                      // union-find
+        ds_repr_(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
@@ -241,7 +247,7 @@ class Persistent_cohomology {
         column_pool_(),  // memory pools for the CAM
         cell_pool_() {
     Simplex_key idx_fil = 0;
-    for (auto & sh : cpx_->filtration_simplex_range()) {
+    for (auto sh : cpx_->filtration_simplex_range()) {
       cpx_->assign_key(sh, idx_fil);
       ++idx_fil;
       dsets_.make_set(cpx_->key(sh));
@@ -264,13 +270,10 @@ class Persistent_cohomology {
   }
 
   ~Persistent_cohomology() {
-// Clean the remaining columns in the matrix.
-    for (auto & cam_ref : cam_) {
-      cam_ref.col_.clear();
-    }
-// Clean the transversal lists
+    // Clean the transversal lists
     for (auto & transverse_ref : transverse_idx_) {
-      transverse_ref.second.row_->clear();
+      // Destruct all the cells
+      transverse_ref.second.row_->clear_and_dispose([&](Cell*p){p->~Cell();});
       delete transverse_ref.second.row_;
     }
   }
@@ -335,18 +338,18 @@ class Persistent_cohomology {
 
       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()));
+        persistent_pairs_.emplace_back(
+            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()));
+      persistent_pairs_.emplace_back(
+          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_));
+      persistent_pairs_.emplace_back(
+          cpx_->simplex(cocycle.first), cpx_->null_simplex(), cocycle.second.characteristics_);
     }
   }
 
@@ -386,8 +389,8 @@ class Persistent_cohomology {
       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()));
+          persistent_pairs_.emplace_back(
+              cpx_->simplex(idx_coc_v), sigma, coeff_field_.characteristic());
         }
         // Maintain the index of the 0-cocycle alive.
         if (kv != idx_coc_v) {
@@ -401,8 +404,8 @@ class Persistent_cohomology {
         }
       } 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()));
+          persistent_pairs_.emplace_back(
+              cpx_->simplex(idx_coc_u), sigma, coeff_field_.characteristic());
         }
         // Maintain the index of the 0-cocycle alive.
         if (ku != idx_coc_u) {
@@ -428,9 +431,12 @@ class Persistent_cohomology {
       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;
+    // with multiplicity. We used to sum the coefficients directly in
+    // annotations_in_boundary by using a map, we now do it later.
+    typedef std::pair<Column *, int> annotation_t;
+    // Danger: not thread-safe!
+    static std::vector<annotation_t> annotations_in_boundary;
+    annotations_in_boundary.clear();
     int sign = 1 - 2 * (dim_sigma % 2);  // \in {-1,1} provides the sign in the
                                          // alternate sum in the boundary.
     Simplex_key key;
@@ -442,22 +448,29 @@ class Persistent_cohomology {
         // 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;
-          }
+	  annotations_in_boundary.emplace_back(curr_col, sign);
         }
       }
       sign = -sign;
     }
+    // Place identical annotations consecutively so we can easily sum their multiplicities.
+    std::sort(annotations_in_boundary.begin(), annotations_in_boundary.end(),
+	[](annotation_t const& a, annotation_t const& b) { return a.first < b.first; });
+
     // 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
+    for (auto ann_it = annotations_in_boundary.begin(); ann_it != annotations_in_boundary.end(); /**/) {
+      Column* col = ann_it->first;
+      int mult = ann_it->second;
+      while (++ann_it != annotations_in_boundary.end() && ann_it->first == col) {
+	mult += ann_it->second;
+      }
+      // The following test is just a heuristic, it is not required, and it is fine that is misses p == 0.
+      if (mult != coeff_field_.additive_identity()) {  // For all columns in the boundary,
+        for (auto cell_ref : col->col_) {  // insert every cell in map_a_ds with multiplicity
+          Arith_element w_y = coeff_field_.times(cell_ref.coefficient_, mult);  // 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));
@@ -525,8 +538,8 @@ class Persistent_cohomology {
                       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));
+    Column * new_col = column_pool_.construct(key);
+    Cell * new_cell = cell_pool_.construct(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
@@ -552,9 +565,9 @@ class Persistent_cohomology {
                        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
+      persistent_pairs_.emplace_back(cpx_->simplex(death_key)  // creator
+          , sigma                                              // destructor
+          , charac);                                           // fields
     }
 
     auto death_key_row = transverse_idx_.find(death_key);  // Find the beginning of the row.
@@ -570,9 +583,8 @@ class Persistent_cohomology {
         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();
+        for (auto& col_cell : curr_col->col_) {
+          col_cell.base_hook_cam_h::unlink();
         }
 
         // Remove the column from the CAM before modifying its value
@@ -587,9 +599,9 @@ class Persistent_cohomology {
           // 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) {
+            for (auto& col_cell : curr_col->col_) {
               // re-establish the row links
-              transverse_idx_[col_cell_it->key_].row_->push_front(*col_cell_it);
+              transverse_idx_[col_cell.key_].row_->push_front(col_cell);
             }
           } else {  // There is already an identical column in the CAM:
             // merge two disjoint sets.
@@ -599,6 +611,8 @@ class Persistent_cohomology {
             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;
+            // intrusive containers don't own their elements, we have to release them manually
+            curr_col->col_.clear_and_dispose([&](Cell*p){cell_pool_.destroy(p);});
             column_pool_.destroy(curr_col);  // delete curr_col;
           }
         }
@@ -694,7 +708,7 @@ class Persistent_cohomology {
   void output_diagram(std::ostream& ostream = std::cout) {
 
     cmp_intervals_by_length cmp(cpx_);
-    persistent_pairs_.sort(cmp);
+    std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
     bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
     for (auto pair : persistent_pairs_) {
       // Special case on windows, inf is "1.#INF" (cf. unitary tests and R package TDA)
@@ -709,13 +723,11 @@ class Persistent_cohomology {
     }
   }
 
-  void write_output_diagram(std::string diagram_name)
-  {
-    std::ofstream           diagram_out(diagram_name.c_str());
-    cmp_intervals_by_length cmp( cpx_ );
-    persistent_pairs_.sort( cmp );
-    for(auto pair : persistent_pairs_)
-    {
+  void write_output_diagram(std::string diagram_name) {
+    std::ofstream diagram_out(diagram_name.c_str());
+    cmp_intervals_by_length cmp(cpx_);
+    std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
+    for (auto pair : persistent_pairs_) {
     diagram_out << cpx_->dimension(get<0>(pair)) << " "
           << cpx_->filtration(get<0>(pair)) << " "
           << cpx_->filtration(get<1>(pair)) << std::endl;
@@ -742,6 +754,7 @@ class Persistent_cohomology {
   Complex_ds * cpx_;
   int dim_max_;
   CoefficientField coeff_field_;
+  size_t num_simplices_;
 
   /*  Disjoint sets data structure to link the model of FilteredComplex
    * with the compressed annotation matrix.
@@ -760,11 +773,11 @@ class Persistent_cohomology {
   /*  Key -> row. */
   std::map<Simplex_key, cocycle> transverse_idx_;
   /* Persistent intervals. */
-  std::list<Persistent_interval> persistent_pairs_;
+  std::vector<Persistent_interval> persistent_pairs_;
   length_interval interval_length_policy;
 
-  boost::object_pool<Column> column_pool_;
-  boost::object_pool<Cell> cell_pool_;
+  Simple_object_pool<Column> column_pool_;
+  Simple_object_pool<Cell> cell_pool_;
 };
 
 /** @} */  // end defgroup persistent_cohomology
@@ -773,4 +786,4 @@ class Persistent_cohomology {
 
 }  // namespace Gudhi
 
-#endif  // SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_H_
+#endif  // 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 2a4c8692..6db16e69 100644
--- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h
+++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h
@@ -20,8 +20,8 @@
  *    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_
+#ifndef PERSISTENT_COHOMOLOGY_FIELD_ZP_H_
+#define PERSISTENT_COHOMOLOGY_FIELD_ZP_H_
 
 #include <utility>
 #include <vector>
@@ -41,9 +41,7 @@ class Field_Zp {
 
   Field_Zp()
       : Prime(0),
-        inverse_(),
-        mult_id_all(1),
-        add_id_all(0) {
+        inverse_() {
   }
 
   void init(int charac) {
@@ -81,14 +79,14 @@ class Field_Zp {
   }
 
   /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/
-  const Element& additive_identity() const {
-    return add_id_all;
+  Element additive_identity() const {
+    return 0;
   }
   /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/
-  const Element& multiplicative_identity(Element = 0) const {
-    return mult_id_all;
+  Element multiplicative_identity(Element = 0) const {
+    return 1;
   }
-  /** Returns the inverse in the field. Modifies P.*/
+  /** 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
@@ -101,7 +99,7 @@ class Field_Zp {
   }
 
   /** \brief Returns the characteristic \f$p\f$ of the field.*/
-  const int& characteristic() const {
+  int characteristic() const {
     return Prime;
   }
 
@@ -109,12 +107,10 @@ class Field_Zp {
   int 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 persistent_cohomology
 
 }  // namespace Gudhi
 
-#endif  // SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_FIELD_ZP_H_
+#endif  // 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 c6fd5282..38bc08d1 100644
--- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h
+++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h
@@ -20,8 +20,8 @@
  *    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_
+#ifndef PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_
+#define PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_
 
 #include <gmpxx.h>
 
@@ -81,8 +81,7 @@ class Multi_field {
     // 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);
+      prod_characteristics_ *= p;
     }
 
     // Uvect_
@@ -142,7 +141,7 @@ class Multi_field {
     return prod_characteristics_;
   }
 
-  /** Returns the inverse in the field. Modifies P.*/
+  /** Returns the inverse in the field. Modifies P. ??? */
   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)
@@ -153,12 +152,12 @@ class Multi_field {
     mpz_invert(inv_qt.get_mpz_t(), x.get_mpz_t(), QT.get_mpz_t());
 
     assert(prod_characteristics_ > 0);  // division by zero + non negative values
-    return std::pair<Element, Element>(
-        (inv_qt * multiplicative_identity(QT)) % prod_characteristics_, QT);
+    return { (inv_qt * multiplicative_identity(QT)) % prod_characteristics_, QT };
   }
   /** Returns -x * y.*/
   Element times_minus(const Element& x, const Element& y) {
     assert(prod_characteristics_ > 0);  // division by zero + non negative values
+    /* This assumes that (x*y)%pc cannot be zero, but Field_Zp has specific code for the 0 case ??? */
     return prod_characteristics_ - ((x * y) % prod_characteristics_);
   }
 
@@ -183,4 +182,4 @@ class Multi_field {
 
 }  // namespace Gudhi
 
-#endif  // SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_
+#endif  // 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 fcec819a..612658e6 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
@@ -20,14 +20,14 @@
  *    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_
+#ifndef PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_
+#define PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_
 
-#include <list>
+#include <boost/tuple/tuple.hpp>
+#include <boost/intrusive/set.hpp>
+#include <boost/intrusive/list.hpp>
 
-#include "boost/tuple/tuple.hpp"
-#include "boost/intrusive/set.hpp"
-#include "boost/intrusive/list.hpp"
+#include <list>
 
 namespace Gudhi {
 
@@ -80,40 +80,33 @@ class Persistent_cohomology_cell : public base_hook_cam_h,
  * 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.
+ *
+ * Movable but not Copyable.
  */
 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;
 
+ public:
   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)
+  explicit Persistent_cohomology_column(SimplexKey key)
       : col_(),
-        class_key_(other.class_key_) {
-    if (!other.col_.empty())
-      std::cerr << "Copying a non-empty column.\n";
-  }
+        class_key_(key) {}
 
   /** \brief Returns true iff the column is null.*/
-  bool is_null() {
+  bool is_null() const {
     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() {
+  SimplexKey class_key() const {
     return class_key_;
   }
 
@@ -145,4 +138,4 @@ class Persistent_cohomology_column : public boost::intrusive::set_base_hook<
 
 }  // namespace Gudhi
 
-#endif  // SRC_PERSISTENT_COHOMOLOGY_INCLUDE_GUDHI_PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_
+#endif  // PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_
diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
index 55bc7066..194b3e74 100644
--- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
+++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp
@@ -1,15 +1,14 @@
-#define BOOST_TEST_MODULE persistent_cohomology 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 <algorithm>
 #include <utility> // std::pair, std::make_pair
-
 #include <cmath> // float comparison
 #include <limits>
 
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "persistent_cohomology"
+#include <boost/test/unit_test.hpp>
+
 #include "gudhi/graph_simplicial_complex.h"
 #include "gudhi/reader_utils.h"
 #include "gudhi/Simplex_tree.h"
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
index 18a4725e..703682e1 100644
--- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp
+++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp
@@ -1,15 +1,14 @@
-#define BOOST_TEST_MODULE persistent_cohomology_multi_field 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 <algorithm>
 #include <utility> // std::pair, std::make_pair
-
 #include <cmath> // float comparison
 #include <limits>
 
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "persistent_cohomology_multi_field"
+#include <boost/test/unit_test.hpp>
+
 #include "gudhi/graph_simplicial_complex.h"
 #include "gudhi/reader_utils.h"
 #include "gudhi/Simplex_tree.h"