Merge last trunk modifications

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

Former-commit-id: 2625074cdd2104ead455753c2f2fb8583be89478

46 files changed, 990 insertions, 166 deletions

diff --git a/src/.clang-format b/src/.clang-format
new file mode 100644
index 00000000..dd4590c8
--- /dev/null
+++ b/src/.clang-format
@@ -0,0 +1,90 @@
+---
+Language:        Cpp
+# BasedOnStyle:  Google
+AccessModifierOffset: -1
+AlignAfterOpenBracket: Align
+AlignConsecutiveAssignments: false
+AlignConsecutiveDeclarations: false
+AlignEscapedNewlinesLeft: true
+AlignOperands:   true
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: false
+AllowShortCaseLabelsOnASingleLine: false
+AllowShortFunctionsOnASingleLine: All
+AllowShortIfStatementsOnASingleLine: true
+AllowShortLoopsOnASingleLine: true
+AlwaysBreakAfterDefinitionReturnType: None
+AlwaysBreakAfterReturnType: None
+AlwaysBreakBeforeMultilineStrings: true
+AlwaysBreakTemplateDeclarations: true
+BinPackArguments: true
+BinPackParameters: true
+BraceWrapping:   
+  AfterClass:      false
+  AfterControlStatement: false
+  AfterEnum:       false
+  AfterFunction:   false
+  AfterNamespace:  false
+  AfterObjCDeclaration: false
+  AfterStruct:     false
+  AfterUnion:      false
+  BeforeCatch:     false
+  BeforeElse:      false
+  IndentBraces:    false
+BreakBeforeBinaryOperators: None
+BreakBeforeBraces: Attach
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+ColumnLimit:     120
+CommentPragmas:  '^ IWYU pragma:'
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
+ConstructorInitializerIndentWidth: 4
+ContinuationIndentWidth: 4
+Cpp11BracedListStyle: true
+DerivePointerAlignment: true
+DisableFormat:   false
+ExperimentalAutoDetectBinPacking: false
+ForEachMacros:   [ foreach, Q_FOREACH, BOOST_FOREACH ]
+IncludeCategories: 
+  - Regex:           '^<.*\.h>'
+    Priority:        1
+  - Regex:           '^<.*'
+    Priority:        2
+  - Regex:           '.*'
+    Priority:        3
+IndentCaseLabels: true
+IndentWidth:     2
+IndentWrappedFunctionNames: false
+KeepEmptyLinesAtTheStartOfBlocks: false
+MacroBlockBegin: ''
+MacroBlockEnd:   ''
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+ObjCBlockIndentWidth: 2
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: false
+PenaltyBreakBeforeFirstCallParameter: 1
+PenaltyBreakComment: 300
+PenaltyBreakFirstLessLess: 120
+PenaltyBreakString: 1000
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 200
+PointerAlignment: Left
+ReflowComments:  true
+SortIncludes:    false
+SpaceAfterCStyleCast: false
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeParens: ControlStatements
+SpaceInEmptyParentheses: false
+SpacesBeforeTrailingComments: 2
+SpacesInAngles:  false
+SpacesInContainerLiterals: true
+SpacesInCStyleCastParentheses: false
+SpacesInParentheses: false
+SpacesInSquareBrackets: false
+Standard:        Auto
+TabWidth:        8
+UseTab:          Never
+...
+
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
index 5a87b9b8..f395de65 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h
@@ -97,7 +97,7 @@ class Bitmap_cubical_complex : public T {
    * with filtration on top dimensional cells.
    **/
   Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
-                         const std::vector<typename T::filtration_type>& top_dimensional_cells) :
+                         const std::vector<Filtration_value>& top_dimensional_cells) :
       T(dimensions, top_dimensional_cells),
       key_associated_to_simplex(this->total_number_of_cells + 1) {
     for (size_t i = 0; i != this->total_number_of_cells; ++i) {
@@ -111,13 +111,13 @@ class Bitmap_cubical_complex : public T {
 
   /**
    * Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells
-   * in the following directions and vector of element of a type T::filtration_type
+   * in the following directions and vector of element of a type Filtration_value
    * with filtration on top dimensional cells. The last parameter of the constructor is a vector of boolean of a length
    * equal to the dimension of cubical complex.
    * If the position i on this vector is true, then we impose periodic boundary conditions in this direction.
    **/
   Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
-                         const std::vector<typename T::filtration_type>& top_dimensional_cells,
+                         const std::vector<Filtration_value>& top_dimensional_cells,
                          std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed) :
       T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed),
       key_associated_to_simplex(this->total_number_of_cells + 1) {
@@ -170,20 +170,20 @@ class Bitmap_cubical_complex : public T {
     if (globalDbg) {
       std::cerr << "unsigned dimension(const Simplex_handle& sh)\n";
     }
-    if (sh != std::numeric_limits<Simplex_handle>::max()) return this->get_dimension_of_a_cell(sh);
+    if (sh != null_simplex()) return this->get_dimension_of_a_cell(sh);
     return -1;
   }
 
   /**
    * Return the filtration of a cell pointed by the Simplex_handle.
    **/
-  typename T::filtration_type filtration(Simplex_handle sh) {
+  Filtration_value filtration(Simplex_handle sh) {
     if (globalDbg) {
-      std::cerr << "T::filtration_type filtration(const Simplex_handle& sh)\n";
+      std::cerr << "Filtration_value filtration(const Simplex_handle& sh)\n";
     }
     // Returns the filtration value of a simplex.
-    if (sh != std::numeric_limits<Simplex_handle>::max()) return this->data[sh];
-    return std::numeric_limits<Simplex_handle>::max();
+    if (sh != null_simplex()) return this->data[sh];
+    return std::numeric_limits<Filtration_value>::infinity();
   }
 
   /**
@@ -203,7 +203,7 @@ class Bitmap_cubical_complex : public T {
     if (globalDbg) {
       std::cerr << "Simplex_key key(const Simplex_handle& sh)\n";
     }
-    if (sh != std::numeric_limits<Simplex_handle>::max()) {
+    if (sh != null_simplex()) {
       return this->key_associated_to_simplex[sh];
     }
     return this->null_key();
@@ -216,7 +216,7 @@ class Bitmap_cubical_complex : public T {
     if (globalDbg) {
       std::cerr << "Simplex_handle simplex(Simplex_key key)\n";
     }
-    if (key != std::numeric_limits<Simplex_handle>::max()) {
+    if (key != null_key()) {
       return this->simplex_associated_to_key[ key ];
     }
     return null_simplex();
@@ -229,7 +229,7 @@ class Bitmap_cubical_complex : public T {
     if (globalDbg) {
       std::cerr << "void assign_key(Simplex_handle& sh, Simplex_key key)\n";
     }
-    if (key == std::numeric_limits<Simplex_handle>::max()) return;
+    if (key == null_key()) return;
     this->key_associated_to_simplex[sh] = key;
     this->simplex_associated_to_key[key] = sh;
   }
@@ -566,8 +566,9 @@ class is_before_in_filtration {
   bool operator()(const typename Bitmap_cubical_complex<T>::Simplex_handle& sh1,
                   const typename Bitmap_cubical_complex<T>::Simplex_handle& sh2) const {
     // Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
-    typename T::filtration_type fil1 = CC_->data[sh1];
-    typename T::filtration_type fil2 = CC_->data[sh2];
+    typedef typename T::filtration_type Filtration_value;
+    Filtration_value fil1 = CC_->data[sh1];
+    Filtration_value fil2 = CC_->data[sh2];
     if (fil1 != fil2) {
       return fil1 < fil2;
     }
diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp
index 238d99ad..1408681a 100644
--- a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp
+++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp
@@ -36,8 +36,8 @@ int main(int argc, char** argv) {
         " distance (set by default to zero). The program will now terminate \n";
     return -1;
   }
-  std::vector<std::pair<double, double>> diag1 = read_persistence_intervals_in_dimension(argv[1]);
-  std::vector<std::pair<double, double>> diag2 = read_persistence_intervals_in_dimension(argv[2]);
+  std::vector<std::pair<double, double>> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]);
+  std::vector<std::pair<double, double>> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]);
 
   double tolerance = 0.;
   if (argc == 4) {
diff --git a/src/Doxyfile b/src/Doxyfile
index 2fecf5fb..6c01aefc 100644
--- a/src/Doxyfile
+++ b/src/Doxyfile
@@ -38,7 +38,7 @@ PROJECT_NAME           = "GUDHI"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = "2.0.0"
+PROJECT_NUMBER         = "2.0.1-rc1"
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
diff --git a/src/Persistent_cohomology/example/rips_distance_matrix_persistence.cpp b/src/Persistent_cohomology/example/rips_distance_matrix_persistence.cpp
index 8517e7f6..d38808c7 100644
--- a/src/Persistent_cohomology/example/rips_distance_matrix_persistence.cpp
+++ b/src/Persistent_cohomology/example/rips_distance_matrix_persistence.cpp
@@ -57,7 +57,7 @@ int main(int argc, char * argv[]) {
 
   program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence);
 
-  Distance_matrix distances = read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file);
+  Distance_matrix distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file);
   Rips_complex rips_complex_from_file(distances, threshold);
 
   // Construct the Rips complex in a Simplex Tree
diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
index 672fda48..e0a147b3 100644
--- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
+++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
@@ -591,10 +591,17 @@ class Persistent_cohomology {
     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);
+    bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
     for (auto pair : persistent_pairs_) {
-    diagram_out << cpx_->dimension(get<0>(pair)) << " "
-          << cpx_->filtration(get<0>(pair)) << " "
-          << cpx_->filtration(get<1>(pair)) << std::endl;
+      // Special case on windows, inf is "1.#INF"
+      if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits<Filtration_value>::infinity()) {
+        diagram_out << cpx_->dimension(get<0>(pair)) << " "
+              << cpx_->filtration(get<0>(pair)) << " inf" << std::endl;
+      } else {
+        diagram_out << cpx_->dimension(get<0>(pair)) << " "
+              << cpx_->filtration(get<0>(pair)) << " "
+              << cpx_->filtration(get<1>(pair)) << std::endl;
+      }
     }
   }
 
diff --git a/src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp b/src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp
index 7ae8126f..9e182f1e 100644
--- a/src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp
+++ b/src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp
@@ -32,7 +32,7 @@ int main(int argc, char **argv) {
   // Init of a Rips complex from a distance matrix in a csv file
   // Default separator is ';'
   // ----------------------------------------------------------------------------
-  Distance_matrix distances = read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_file_name);
+  Distance_matrix distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_file_name);
   Rips_complex rips_complex_from_file(distances, threshold);
 
   std::streambuf* streambufffer;
diff --git a/src/Rips_complex/test/test_rips_complex.cpp b/src/Rips_complex/test/test_rips_complex.cpp
index fc2179f2..fc83f5f7 100644
--- a/src/Rips_complex/test/test_rips_complex.cpp
+++ b/src/Rips_complex/test/test_rips_complex.cpp
@@ -244,7 +244,7 @@ BOOST_AUTO_TEST_CASE(Rips_doc_csv_file) {
   std::cout << "========== CSV FILE NAME = " << csv_file_name << " - Rips threshold=" <<
       rips_threshold << "==========" << std::endl;
 
-  Distance_matrix distances = read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_file_name);
+  Distance_matrix distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_file_name);
   Rips_complex rips_complex_from_file(distances, rips_threshold);
 
   const int DIMENSION_1 = 1;
diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake
index 8f486118..dbf2106a 100644
--- a/src/cmake/modules/GUDHI_third_party_libraries.cmake
+++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake
@@ -117,15 +117,13 @@ if(NOT GUDHI_CYTHON_PATH)
 endif(NOT GUDHI_CYTHON_PATH)
 
 if(PYTHONINTERP_FOUND AND CYTHON_FOUND)
-  # Unitary tests are available through py.test
-  find_program( PYTEST_PATH py.test )
   # Default found version 2
   if(PYTHON_VERSION_MAJOR EQUAL 2)
     # Documentation generation is available through sphinx
     find_program( SPHINX_PATH sphinx-build )
   elseif(PYTHON_VERSION_MAJOR EQUAL 3)
     # No sphinx-build in Pyton3, just hack it
-    set(SPHINX_PATH "${CMAKE_SOURCE_DIR}/${GUDHI_CYTHON_PATH}/doc/python3-sphinx-build")
+    set(SPHINX_PATH "${CMAKE_CURRENT_BINARY_DIR}/${GUDHI_CYTHON_PATH}/doc/python3-sphinx-build")
   else()
     message(FATAL_ERROR "ERROR: Try to compile the Cython interface. Python version ${PYTHON_VERSION_STRING} is not valid.")
   endif(PYTHON_VERSION_MAJOR EQUAL 2)
diff --git a/src/common/doc/file_formats.h b/src/common/doc/file_formats.h
index c145b271..d715aa4d 100644
--- a/src/common/doc/file_formats.h
+++ b/src/common/doc/file_formats.h
@@ -2,7 +2,7 @@
 *    (Geometric Understanding in Higher Dimensions) is a generic C++
 *    library for computational topology.
 *
-*    Author(s):       Clément Jamin
+*    Author(s):       Clément Jamin
 *
 *    Copyright (C) 2017  INRIA
 *
@@ -26,7 +26,7 @@
 namespace Gudhi {
 
 /*! \page fileformats File formats
- 
+
  \tableofcontents
 
  \section FileFormatsPers Persistence Diagram
@@ -34,20 +34,25 @@ namespace Gudhi {
  Such a file, whose extension is usually `.pers`, contains a list of persistence intervals.<br>
  Lines starting with `#` are ignored (comments).<br>
  Other lines might contain 2, 3 or 4 values (the number of values on each line must be the same for all lines):
- \code{.unparsed}
+ \verbatim
    [[field] dimension] birth death
- \endcode
+ \endverbatim
 
  Here is a simple sample file:
- \code{.unparsed}
-   # Beautiful persistence diagram
+ \verbatim
+   # Persistence diagram example
    2 2.7 3.7
    2 9.6 14.
+   # Some comments
    3 34.2 34.974
    4 3. inf
- \endcode
+ \endverbatim
 
  Other sample files can be found in the `data/persistence_diagram` folder.
+
+ Such files can be generated with `Gudhi::persistent_cohomology::Persistent_cohomology::output_diagram()` and read with
+ `Gudhi::read_persistence_intervals_and_dimension()`, `Gudhi::read_persistence_intervals_grouped_by_dimension()` or
+ `Gudhi::read_persistence_intervals_in_dimension()`.
 */
 }  // namespace Gudhi
 
diff --git a/src/common/include/gudhi/reader_utils.h b/src/common/include/gudhi/reader_utils.h
index f1684d78..90be4fc7 100644
--- a/src/common/include/gudhi/reader_utils.h
+++ b/src/common/include/gudhi/reader_utils.h
@@ -1,5 +1,5 @@
-/*    This file is part of the Gudhi Library. The Gudhi library 
- *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
+/*    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):       Clement Maria, Pawel Dlotko, Clement Jamin
@@ -36,6 +36,9 @@
 #include <string>
 #include <vector>
 #include <utility>  // for pair
+#include <tuple>  // for std::make_tuple
+
+namespace Gudhi {
 
 // Keep this file tag for Doxygen to parse the code, otherwise, functions are not documented.
 // It is required for global functions and variables.
@@ -52,7 +55,7 @@
  * X21 X22 ... X2d<br>
  * etc<br>
  */
-inline void read_points(std::string file_name, std::vector< std::vector< double > > & points) {
+inline void read_points(std::string file_name, std::vector<std::vector<double>>& points) {
   std::ifstream in_file(file_name.c_str(), std::ios::in);
   if (!in_file.is_open()) {
     std::cerr << "Unable to open file " << file_name << std::endl;
@@ -62,14 +65,13 @@ inline void read_points(std::string file_name, std::vector< std::vector< double
   std::string line;
   double x;
   while (getline(in_file, line)) {
-    std::vector< double > point;
+    std::vector<double> point;
     std::istringstream iss(line);
     while (iss >> x) {
       point.push_back(x);
     }
     // Check for empty lines
-    if (!point.empty())
-      points.push_back(point);
+    if (!point.empty()) points.push_back(point);
   }
   in_file.close();
 }
@@ -90,7 +92,7 @@ inline void read_points(std::string file_name, std::vector< std::vector< double
  * Every simplex must appear exactly once.
  * Simplices of dimension more than 1 are ignored.
  */
-template< typename Graph_t, typename Filtration_value, typename Vertex_handle >
+template <typename Graph_t, typename Filtration_value, typename Vertex_handle>
 Graph_t read_graph(std::string file_name) {
   std::ifstream in_(file_name.c_str(), std::ios::in);
   if (!in_.is_open()) {
@@ -100,10 +102,10 @@ Graph_t read_graph(std::string file_name) {
     throw std::invalid_argument(error_str);
   }
 
-  typedef std::pair< Vertex_handle, Vertex_handle > Edge_t;
-  std::vector< Edge_t > edges;
-  std::vector< Filtration_value > edges_fil;
-  std::map< Vertex_handle, Filtration_value > vertices;
+  typedef std::pair<Vertex_handle, Vertex_handle> Edge_t;
+  std::vector<Edge_t> edges;
+  std::vector<Filtration_value> edges_fil;
+  std::map<Vertex_handle, Filtration_value> vertices;
 
   std::string line;
   int dim;
@@ -113,8 +115,7 @@ Graph_t read_graph(std::string file_name) {
     std::istringstream iss(line);
     while (iss >> dim) {
       switch (dim) {
-        case 0:
-        {
+        case 0: {
           iss >> u;
           iss >> fil;
           vertices[u] = fil;
@@ -123,8 +124,7 @@ Graph_t read_graph(std::string file_name) {
           }
           break;
         }
-        case 1:
-        {
+        case 1: {
           iss >> u;
           iss >> v;
           iss >> fil;
@@ -132,16 +132,13 @@ Graph_t read_graph(std::string file_name) {
           edges_fil.push_back(fil);
           break;
         }
-        default:
-        {
-          break;
-        }
+        default: { break; }
       }
     }
   }
   in_.close();
 
-  if ((size_t) (max_h + 1) != vertices.size()) {
+  if ((size_t)(max_h + 1) != vertices.size()) {
     std::cerr << "Error: vertices must be labeled from 0 to n-1 \n";
   }
 
@@ -169,8 +166,8 @@ Graph_t read_graph(std::string file_name) {
  * Every simplex must appear exactly once.
  * Simplices of dimension more than 1 are ignored.
  */
-template< typename Vertex_handle, typename Filtration_value >
-bool read_simplex(std::istream & in_, std::vector< Vertex_handle > & simplex, Filtration_value & fil) {
+template <typename Vertex_handle, typename Filtration_value>
+bool read_simplex(std::istream& in_, std::vector<Vertex_handle>& simplex, Filtration_value& fil) {
   int dim = 0;
   if (!(in_ >> dim)) return false;
   Vertex_handle v;
@@ -194,8 +191,8 @@ bool read_simplex(std::istream & in_, std::vector< Vertex_handle > & simplex, Fi
  * The key of a simplex is its position in the filtration order and also the number of its row in the file.
  * Dimi ki1 ki2 ... kiDimi Fili means that the ith simplex in the filtration has dimension Dimi, filtration value
  * fil1 and simplices with key ki1 ... kiDimi in its boundary.*/
-template< typename Simplex_key, typename Filtration_value >
-bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundary, Filtration_value & fil) {
+template <typename Simplex_key, typename Filtration_value>
+bool read_hasse_simplex(std::istream& in_, std::vector<Simplex_key>& boundary, Filtration_value& fil) {
   int dim;
   if (!(in_ >> dim)) return false;
   if (dim == 0) {
@@ -214,7 +211,7 @@ bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundar
 /**
  * @brief Read a lower triangular distance matrix from a csv file. We assume that the .csv store the whole
  * (square) matrix.
- * 
+ *
  * @author Pawel Dlotko
  *
  * Square matrix file format:<br>
@@ -231,13 +228,13 @@ bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundar
  * Dj1;Dj2;...;Dj(j-1);<br>
  *
  **/
-template< typename Filtration_value >
-std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from_csv_file(const std::string& filename,
-                                                                                          const char separator = ';') {
+template <typename Filtration_value>
+std::vector<std::vector<Filtration_value>> read_lower_triangular_matrix_from_csv_file(const std::string& filename,
+                                                                                      const char separator = ';') {
 #ifdef DEBUG_TRACES
   std::cout << "Using procedure read_lower_triangular_matrix_from_csv_file \n";
 #endif  // DEBUG_TRACES
-  std::vector< std::vector< Filtration_value > > result;
+  std::vector<std::vector<Filtration_value>> result;
   std::ifstream in;
   in.open(filename.c_str());
   if (!in.is_open()) {
@@ -248,7 +245,7 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
 
   // the first line is emtpy, so we ignore it:
   std::getline(in, line);
-  std::vector< Filtration_value > values_in_this_line;
+  std::vector<Filtration_value> values_in_this_line;
   result.push_back(values_in_this_line);
 
   int number_of_line = 0;
@@ -256,11 +253,10 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
   // first, read the file line by line to a string:
   while (std::getline(in, line)) {
     // if line is empty, break
-    if (line.size() == 0)
-      break;
+    if (line.size() == 0) break;
 
     // if the last element of a string is comma:
-    if (line[ line.size() - 1 ] == separator) {
+    if (line[line.size() - 1] == separator) {
       // then shrink the string by one
       line.pop_back();
     }
@@ -273,7 +269,7 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
     // and now read the doubles.
 
     int number_of_entry = 0;
-    std::vector< Filtration_value > values_in_this_line;
+    std::vector<Filtration_value> values_in_this_line;
     while (iss.good()) {
       double entry;
       iss >> entry;
@@ -282,7 +278,7 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
       }
       ++number_of_entry;
     }
-    if (!values_in_this_line.empty())result.push_back(values_in_this_line);
+    if (!values_in_this_line.empty()) result.push_back(values_in_this_line);
     ++number_of_line;
   }
   in.close();
@@ -309,7 +305,6 @@ Note: the function does not check that birth <= death.
 **/
 template <typename OutputIterator>
 void read_persistence_intervals_and_dimension(std::string const& filename, OutputIterator out) {
-
   std::ifstream in(filename);
   if (!in.is_open()) {
     std::string error_str("read_persistence_intervals_and_dimension - Unable to open file ");
@@ -325,13 +320,12 @@ void read_persistence_intervals_and_dimension(std::string const& filename, Outpu
       double numbers[4];
       int n = sscanf(line.c_str(), "%lf %lf %lf %lf", &numbers[0], &numbers[1], &numbers[2], &numbers[3]);
       if (n >= 2) {
-        //int field = (n == 4 ? static_cast<int>(numbers[0]) : -1);
         int dim = (n >= 3 ? static_cast<int>(numbers[n - 3]) : -1);
         *out++ = std::make_tuple(dim, numbers[n - 2], numbers[n - 1]);
       }
     }
   }
-} // read_persistence_diagram_from_file
+}
 
 /**
 Reads a file containing persistence intervals.
@@ -340,33 +334,36 @@ The return value is an `std::map<dim, std::vector<std::pair<birth, death>>>`
 where `dim` is an `int`, `birth` a `double`, and `death` a `double`.
 Note: the function does not check that birth <= death.
 **/
-inline std::map<int, std::vector<std::pair<double, double>>> read_persistence_intervals_grouped_by_dimension(std::string const& filename) {
-
+inline std::map<int, std::vector<std::pair<double, double>>> read_persistence_intervals_grouped_by_dimension(
+    std::string const& filename) {
   std::map<int, std::vector<std::pair<double, double>>> ret;
   read_persistence_intervals_and_dimension(
-    filename,
-    boost::make_function_output_iterator([&ret](std::tuple<int, double, double> t) { ret[get<0>(t)].push_back(std::make_pair(get<1>(t), get<2>(t))); }));
+      filename, boost::make_function_output_iterator([&ret](std::tuple<int, double, double> t) {
+        ret[get<0>(t)].push_back(std::make_pair(get<1>(t), get<2>(t)));
+      }));
   return ret;
-} // read_persistence_diagram_from_file
-
+}
 
 /**
 Reads a file containing persistence intervals.
 Each line might contain 2, 3 or 4 values: [[field] dimension] birth death
 If `only_this_dim` = -1, dimension is ignored and all lines are returned.
-If `only_this_dim` is >= 0, only the lines where dimension = `only_this_dim` 
+If `only_this_dim` is >= 0, only the lines where dimension = `only_this_dim`
 (or where dimension is not specified) are returned.
 The return value is an `std::vector<std::pair<birth, death>>`
 where `dim` is an `int`, `birth` a `double`, and `death` a `double`.
 Note: the function does not check that birth <= death.
 **/
-inline std::vector<std::pair<double, double>> read_persistence_intervals_in_dimension(std::string const& filename, int only_this_dim = -1) {
-
+inline std::vector<std::pair<double, double>> read_persistence_intervals_in_dimension(std::string const& filename,
+                                                                                      int only_this_dim = -1) {
   std::vector<std::pair<double, double>> ret;
   read_persistence_intervals_and_dimension(
-    filename, 
-    boost::make_function_output_iterator([&ret](std::tuple<int, double, double> t) { ret.emplace_back(get<1>(t), get<2>(t)); }));
+      filename, boost::make_function_output_iterator([only_this_dim, &ret](std::tuple<int, double, double> t) {
+        if (only_this_dim == get<0>(t) || only_this_dim == -1) ret.emplace_back(get<1>(t), get<2>(t));
+      }));
   return ret;
-} // read_persistence_diagram_from_file
+}
+
+}  // namespace Gudhi
 
 #endif  // READER_UTILS_H_
diff --git a/src/common/test/CMakeLists.txt b/src/common/test/CMakeLists.txt
index 5aa426d7..de3e765a 100644
--- a/src/common/test/CMakeLists.txt
+++ b/src/common/test/CMakeLists.txt
@@ -9,10 +9,17 @@ target_link_libraries(Common_test_points_off_reader ${Boost_UNIT_TEST_FRAMEWORK_
 add_executable ( Common_test_distance_matrix_reader test_distance_matrix_reader.cpp )
 target_link_libraries(Common_test_distance_matrix_reader ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
 
+add_executable ( Common_test_persistence_intervals_reader test_persistence_intervals_reader.cpp )
+target_link_libraries(Common_test_persistence_intervals_reader ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+
 # Do not forget to copy test files in current binary dir
 file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
 file(COPY "${CMAKE_SOURCE_DIR}/data/distance_matrix/lower_triangular_distance_matrix.csv" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
 file(COPY "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+file(COPY "${CMAKE_SOURCE_DIR}/src/common/test/persistence_intervals_with_dimension.pers" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+file(COPY "${CMAKE_SOURCE_DIR}/src/common/test/persistence_intervals_with_field.pers" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+file(COPY "${CMAKE_SOURCE_DIR}/src/common/test/persistence_intervals_without_dimension.pers" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
 
 gudhi_add_coverage_test(Common_test_points_off_reader)
 gudhi_add_coverage_test(Common_test_distance_matrix_reader)
+gudhi_add_coverage_test(Common_test_persistence_intervals_reader)
diff --git a/src/common/test/README b/src/common/test/README
index f2a7eb5a..a8e6efe9 100644
--- a/src/common/test/README
+++ b/src/common/test/README
@@ -7,7 +7,7 @@ make
 To launch with details:
 ***********************
 
-./dtoffrw_UT --report_level=detailed --log_level=all 
+./Common_test_points_off_reader --report_level=detailed --log_level=all
 
    ==> echo $? returns 0 in case of success (non-zero otherwise)
 
diff --git a/src/common/test/persistence_intervals_with_dimension.pers b/src/common/test/persistence_intervals_with_dimension.pers
new file mode 100644
index 00000000..406748c8
--- /dev/null
+++ b/src/common/test/persistence_intervals_with_dimension.pers
@@ -0,0 +1,5 @@
+# Simple persistence diagram with dimension
+0 2.7 3.7
+1 9.6 14.
+3 34.2 34.974
+1 3. inf
diff --git a/src/common/test/persistence_intervals_with_field.pers b/src/common/test/persistence_intervals_with_field.pers
new file mode 100644
index 00000000..41dd9f1e
--- /dev/null
+++ b/src/common/test/persistence_intervals_with_field.pers
@@ -0,0 +1,4 @@
+3 0 2.7 3.7
+3 1 9.6 14.
+3 3 34.2 34.974
+3 1 3. inf
diff --git a/src/common/test/persistence_intervals_without_dimension.pers b/src/common/test/persistence_intervals_without_dimension.pers
new file mode 100644
index 00000000..76fa27f3
--- /dev/null
+++ b/src/common/test/persistence_intervals_without_dimension.pers
@@ -0,0 +1,7 @@
+# Simple persistence diagram without dimension
+2.7 3.7
+9.6 14.
+# Another comment
+34.2 34.974
+3. inf
+# End of file
diff --git a/src/common/test/test_distance_matrix_reader.cpp b/src/common/test/test_distance_matrix_reader.cpp
index 95a73bd9..656e6f2e 100644
--- a/src/common/test/test_distance_matrix_reader.cpp
+++ b/src/common/test/test_distance_matrix_reader.cpp
@@ -36,7 +36,7 @@ BOOST_AUTO_TEST_CASE( lower_triangular_distance_matrix )
 {
   Distance_matrix from_lower_triangular;
   // Read lower_triangular_distance_matrix.csv file where the separator is a ','
-  from_lower_triangular = read_lower_triangular_matrix_from_csv_file<double>("lower_triangular_distance_matrix.csv",
+  from_lower_triangular = Gudhi::read_lower_triangular_matrix_from_csv_file<double>("lower_triangular_distance_matrix.csv",
                                                                              ',');
   for (auto& i : from_lower_triangular) {
     for (auto j : i) {
@@ -69,7 +69,7 @@ BOOST_AUTO_TEST_CASE( full_square_distance_matrix )
 {
   Distance_matrix from_full_square;
   // Read full_square_distance_matrix.csv file where the separator is the default one ';'
-  from_full_square = read_lower_triangular_matrix_from_csv_file<double>("full_square_distance_matrix.csv");
+  from_full_square = Gudhi::read_lower_triangular_matrix_from_csv_file<double>("full_square_distance_matrix.csv");
   for (auto& i : from_full_square) {
     for (auto j : i) {
       std::cout << j << " ";
diff --git a/src/common/test/test_persistence_intervals_reader.cpp b/src/common/test/test_persistence_intervals_reader.cpp
new file mode 100644
index 00000000..be299376
--- /dev/null
+++ b/src/common/test/test_persistence_intervals_reader.cpp
@@ -0,0 +1,322 @@
+/*    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) 2017  INRIA
+ *
+ *    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 <iostream>
+#include <vector>
+#include <utility>  // for pair
+#include <tuple>
+#include <limits>  // for inf
+#include <map>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "persistence_intervals_reader"
+#include <boost/test/unit_test.hpp>
+
+using Persistence_intervals_by_dimension = std::map<int, std::vector<std::pair<double, double>>>;
+using Persistence_intervals = std::vector<std::pair<double, double>>;
+// Test files with only 2 parameters (persistence birth and death) per line in file
+BOOST_AUTO_TEST_CASE( persistence_intervals_without_dimension )
+{
+  Persistence_intervals_by_dimension expected_intervals_by_dimension;
+  expected_intervals_by_dimension[-1].push_back(std::make_pair(2.7, 3.7));
+  expected_intervals_by_dimension[-1].push_back(std::make_pair(9.6, 14.));
+  expected_intervals_by_dimension[-1].push_back(std::make_pair(34.2, 34.974));
+  expected_intervals_by_dimension[-1].push_back(std::make_pair(3., std::numeric_limits<double>::infinity()));
+
+  Persistence_intervals_by_dimension persistence_intervals_by_dimension =
+      Gudhi::read_persistence_intervals_grouped_by_dimension("persistence_intervals_without_dimension.pers");
+
+  std::cout << "\nread_persistence_intervals_grouped_by_dimension - expected\n";
+  for (auto map_iter : expected_intervals_by_dimension) {
+    std::cout << "key=" << map_iter.first;
+    for (auto vec_iter : map_iter.second)
+      std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+  }
+
+  std::cout << "\nread_persistence_intervals_grouped_by_dimension - read\n";
+  for (auto map_iter : persistence_intervals_by_dimension) {
+    std::cout << "key=" << map_iter.first;
+    for (auto vec_iter : map_iter.second)
+      std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+  }
+
+  BOOST_CHECK(persistence_intervals_by_dimension == expected_intervals_by_dimension);
+
+  Persistence_intervals expected_intervals_in_dimension;
+  expected_intervals_in_dimension.push_back(std::make_pair(2.7, 3.7));
+  expected_intervals_in_dimension.push_back(std::make_pair(9.6, 14.));
+  expected_intervals_in_dimension.push_back(std::make_pair(34.2, 34.974));
+  expected_intervals_in_dimension.push_back(std::make_pair(3., std::numeric_limits<double>::infinity()));
+
+  Persistence_intervals persistence_intervals_in_dimension =
+      Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_without_dimension.pers");
+
+  std::cout << "\nread_persistence_intervals_in_dimension - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension - read\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_without_dimension.pers", 0);
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_without_dimension.pers", 1);
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_without_dimension.pers", 2);
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_without_dimension.pers", 3);
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+}
+// Test files with 3 parameters (dimension birth death) per line in file
+BOOST_AUTO_TEST_CASE( persistence_intervals_with_dimension )
+{
+  Persistence_intervals_by_dimension expected_intervals_by_dimension;
+  expected_intervals_by_dimension[0].push_back(std::make_pair(2.7, 3.7));
+  expected_intervals_by_dimension[1].push_back(std::make_pair(9.6, 14.));
+  expected_intervals_by_dimension[3].push_back(std::make_pair(34.2, 34.974));
+  expected_intervals_by_dimension[1].push_back(std::make_pair(3., std::numeric_limits<double>::infinity()));
+
+  Persistence_intervals_by_dimension persistence_intervals_by_dimension =
+      Gudhi::read_persistence_intervals_grouped_by_dimension("persistence_intervals_with_dimension.pers");
+
+  std::cout << "\nread_persistence_intervals_grouped_by_dimension - expected\n";
+  for (auto map_iter : expected_intervals_by_dimension) {
+    std::cout << "key=" << map_iter.first;
+    for (auto vec_iter : map_iter.second)
+      std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+  }
+
+  std::cout << "\nread_persistence_intervals_grouped_by_dimension - read\n";
+  for (auto map_iter : persistence_intervals_by_dimension) {
+    std::cout << "key=" << map_iter.first;
+    for (auto vec_iter : map_iter.second)
+      std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+  }
+
+  BOOST_CHECK(persistence_intervals_by_dimension == expected_intervals_by_dimension);
+
+  Persistence_intervals expected_intervals_in_dimension;
+  expected_intervals_in_dimension.push_back(std::make_pair(2.7, 3.7));
+  expected_intervals_in_dimension.push_back(std::make_pair(9.6, 14.));
+  expected_intervals_in_dimension.push_back(std::make_pair(34.2, 34.974));
+  expected_intervals_in_dimension.push_back(std::make_pair(3., std::numeric_limits<double>::infinity()));
+
+  Persistence_intervals persistence_intervals_in_dimension =
+      Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_dimension.pers");
+
+  std::cout << "\nread_persistence_intervals_in_dimension - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  expected_intervals_in_dimension.push_back(std::make_pair(2.7, 3.7));
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_dimension.pers", 0);
+
+  std::cout << "\nread_persistence_intervals_in_dimension 0 - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension 0 - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  expected_intervals_in_dimension.push_back(std::make_pair(9.6, 14.));
+  expected_intervals_in_dimension.push_back(std::make_pair(3., std::numeric_limits<double>::infinity()));
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_dimension.pers", 1);
+
+  std::cout << "\nread_persistence_intervals_in_dimension 1 - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension 1 - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_dimension.pers", 2);
+
+  std::cout << "\nread_persistence_intervals_in_dimension 2 - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension 2 - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  expected_intervals_in_dimension.push_back(std::make_pair(34.2, 34.974));
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_dimension.pers", 3);
+
+  std::cout << "\nread_persistence_intervals_in_dimension 3 - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension 3 - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+}
+
+// Test files with 4 parameters (field dimension birth death) per line in file
+BOOST_AUTO_TEST_CASE( persistence_intervals_with_field )
+{
+  Persistence_intervals_by_dimension expected_intervals_by_dimension;
+  expected_intervals_by_dimension[0].push_back(std::make_pair(2.7, 3.7));
+  expected_intervals_by_dimension[1].push_back(std::make_pair(9.6, 14.));
+  expected_intervals_by_dimension[3].push_back(std::make_pair(34.2, 34.974));
+  expected_intervals_by_dimension[1].push_back(std::make_pair(3., std::numeric_limits<double>::infinity()));
+
+  Persistence_intervals_by_dimension persistence_intervals_by_dimension =
+      Gudhi::read_persistence_intervals_grouped_by_dimension("persistence_intervals_with_field.pers");
+
+  std::cout << "\nread_persistence_intervals_grouped_by_dimension - expected\n";
+  for (auto map_iter : expected_intervals_by_dimension) {
+    std::cout << "key=" << map_iter.first;
+    for (auto vec_iter : map_iter.second)
+      std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+  }
+
+  std::cout << "\nread_persistence_intervals_grouped_by_dimension - read\n";
+  for (auto map_iter : persistence_intervals_by_dimension) {
+    std::cout << "key=" << map_iter.first;
+    for (auto vec_iter : map_iter.second)
+      std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+  }
+
+  BOOST_CHECK(persistence_intervals_by_dimension == expected_intervals_by_dimension);
+
+  Persistence_intervals expected_intervals_in_dimension;
+  expected_intervals_in_dimension.push_back(std::make_pair(2.7, 3.7));
+  expected_intervals_in_dimension.push_back(std::make_pair(9.6, 14.));
+  expected_intervals_in_dimension.push_back(std::make_pair(34.2, 34.974));
+  expected_intervals_in_dimension.push_back(std::make_pair(3., std::numeric_limits<double>::infinity()));
+
+  Persistence_intervals persistence_intervals_in_dimension =
+      Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_field.pers");
+
+  std::cout << "\nread_persistence_intervals_in_dimension - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  expected_intervals_in_dimension.push_back(std::make_pair(2.7, 3.7));
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_field.pers", 0);
+
+  std::cout << "\nread_persistence_intervals_in_dimension 0 - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension 0 - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  expected_intervals_in_dimension.push_back(std::make_pair(9.6, 14.));
+  expected_intervals_in_dimension.push_back(std::make_pair(3., std::numeric_limits<double>::infinity()));
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_field.pers", 1);
+
+  std::cout << "\nread_persistence_intervals_in_dimension 1 - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension 1 - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_field.pers", 2);
+
+  std::cout << "\nread_persistence_intervals_in_dimension 2 - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension 2 - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+  expected_intervals_in_dimension.clear();
+  expected_intervals_in_dimension.push_back(std::make_pair(34.2, 34.974));
+  persistence_intervals_in_dimension =
+        Gudhi::read_persistence_intervals_in_dimension("persistence_intervals_with_field.pers", 3);
+
+  std::cout << "\nread_persistence_intervals_in_dimension 3 - expected\n";
+  for (auto vec_iter : expected_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  std::cout << "\nread_persistence_intervals_in_dimension 3 - read\n";
+  for (auto vec_iter : persistence_intervals_in_dimension)
+    std::cout << " [" << vec_iter.first << " ," << vec_iter.second << "] ";
+
+  BOOST_CHECK(persistence_intervals_in_dimension == expected_intervals_in_dimension);
+
+}
diff --git a/src/cython/CMakeLists.txt b/src/cython/CMakeLists.txt
index 99badffb..ab5ea9ef 100644
--- a/src/cython/CMakeLists.txt
+++ b/src/cython/CMakeLists.txt
@@ -16,7 +16,7 @@ function( add_gudhi_cython_lib THE_LIB )
 endfunction( add_gudhi_cython_lib )
 
 if(CYTHON_FOUND)
-  message("++ ${PYTHON_EXECUTABLE} v.${PYTHON_VERSION_STRING} - Cython is ${CYTHON_EXECUTABLE} - py.test is ${PYTEST_PATH} - Sphinx is ${SPHINX_PATH}")
+  message("++ ${PYTHON_EXECUTABLE} v.${PYTHON_VERSION_STRING} - Cython is ${CYTHON_EXECUTABLE} - Sphinx is ${SPHINX_PATH}")
 
   set(GUDHI_CYTHON_EXTRA_COMPILE_ARGS "${GUDHI_CYTHON_EXTRA_COMPILE_ARGS}'-DBOOST_RESULT_OF_USE_DECLTYPE', ")
   set(GUDHI_CYTHON_EXTRA_COMPILE_ARGS "${GUDHI_CYTHON_EXTRA_COMPILE_ARGS}'-DBOOST_ALL_NO_LIB', ")
@@ -100,7 +100,7 @@ if(CYTHON_FOUND)
   file(COPY "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/doc/")
   # Persistence graphical tools examples
   file(COPY "${CMAKE_SOURCE_DIR}/data/bitmap/3d_torus.txt" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/doc/")
-  file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/doc/")
+  file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/doc/")
 
   if (NOT CGAL_VERSION VERSION_LESS 4.8.1)
     # If CGAL_VERSION >= 4.8.1, include subsampling
@@ -194,6 +194,7 @@ if(CYTHON_FOUND)
   endif(UNIX)
 
   # set sphinx-build in make files
+  configure_file(doc/python3-sphinx-build.in "${CMAKE_CURRENT_BINARY_DIR}/doc/python3-sphinx-build" @ONLY)
   configure_file(doc/Makefile.in "${CMAKE_CURRENT_BINARY_DIR}/doc/Makefile" @ONLY)
   configure_file(doc/make.bat.in "${CMAKE_CURRENT_BINARY_DIR}/doc/make.bat" @ONLY)
 
@@ -305,13 +306,11 @@ if(CYTHON_FOUND)
   set_tests_properties(witness_complex_from_nearest_landmark_table_py_test PROPERTIES ENVIRONMENT "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}")
 
   # Unitary tests are available through py.test
-  if(PYTEST_PATH)
-    add_test(
-        NAME gudhi_cython_py_test
-        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-        COMMAND ${PYTHON_EXECUTABLE} "${PYTEST_PATH}")
-    set_tests_properties(gudhi_cython_py_test PROPERTIES ENVIRONMENT "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}")
-  endif(PYTEST_PATH)
+  add_test(
+      NAME gudhi_cython_py_test
+      WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+      COMMAND ${PYTHON_EXECUTABLE} -m pytest)
+  set_tests_properties(gudhi_cython_py_test PROPERTIES ENVIRONMENT "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}")
 
   # Documentation generation is available through sphinx
   if(SPHINX_PATH)
diff --git a/src/cython/cython/persistence_graphical_tools.py b/src/cython/cython/persistence_graphical_tools.py
index a984633e..fb837e29 100755
--- a/src/cython/cython/persistence_graphical_tools.py
+++ b/src/cython/cython/persistence_graphical_tools.py
@@ -1,11 +1,12 @@
 import matplotlib.pyplot as plt
 import numpy as np
+import os
 
 """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
+   Author(s):       Vincent Rouvreau, Bertrand Michel
 
    Copyright (C) 2016 INRIA
 
@@ -23,15 +24,17 @@ import numpy as np
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
-__author__ = "Vincent Rouvreau"
+__author__ = "Vincent Rouvreau, Bertrand Michel"
 __copyright__ = "Copyright (C) 2016 INRIA"
 __license__ = "GPL v3"
 
-def __min_birth_max_death(persistence):
+def __min_birth_max_death(persistence, band_boot=0.):
     """This function returns (min_birth, max_death) from the persistence.
 
     :param persistence: The persistence to plot.
     :type persistence: list of tuples(dimension, tuple(birth, death)).
+    :param band_boot: bootstrap band
+    :type band_boot: float.
     :returns: (float, float) -- (min_birth, max_death).
     """
     # Look for minimum birth date and maximum death date for plot optimisation
@@ -45,6 +48,8 @@ def __min_birth_max_death(persistence):
             max_death = float(interval[1][0])
         if float(interval[1][0]) < min_birth:
             min_birth = float(interval[1][0])
+    if band_boot > 0.:
+        max_death += band_boot
     return (min_birth, max_death)
 
 """
@@ -59,7 +64,7 @@ def show_palette_values(alpha=0.6):
 
     :param alpha: alpha value in [0.0, 1.0] for horizontal bars (default is 0.6).
     :type alpha: float.
-    :returns: plot -- An horizontal bar plot of dimensions color.
+    :returns: plot the dimension palette values.
     """
     colors = []
     for color in palette:
@@ -70,18 +75,38 @@ def show_palette_values(alpha=0.6):
     plt.barh(y_pos, y_pos + 1, align='center', alpha=alpha, color=colors)
     plt.ylabel('Dimension')
     plt.title('Dimension palette values')
+    return plt
 
-    plt.show()
-
-def plot_persistence_barcode(persistence, alpha=0.6):
+def plot_persistence_barcode(persistence=[], persistence_file='', alpha=0.6, max_barcodes=0):
     """This function plots the persistence bar code.
 
     :param persistence: The persistence to plot.
     :type persistence: list of tuples(dimension, tuple(birth, death)).
+    :param persistence_file: A persistence file style name (reset persistence if both are set).
+    :type persistence_file: string
     :param alpha: alpha value in [0.0, 1.0] for horizontal bars (default is 0.6).
     :type alpha: float.
+    :param max_barcodes: number of maximal barcodes to be displayed
+        (persistence will be sorted by life time if max_barcodes is set)
+    :type max_barcodes: int.
     :returns: plot -- An horizontal bar plot of persistence.
     """
+    if persistence_file is not '':
+        if os.path.isfile(persistence_file):
+            # Reset persistence
+            persistence = []
+            diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file)
+            for key in diag.keys():
+                for persistence_interval in diag[key]:
+                    persistence.append((key, persistence_interval))
+        else:
+            print("file " + persistence_file + " not found.")
+            return None
+
+    if max_barcodes > 0 and max_barcodes < len(persistence):
+        # Sort by life time, then takes only the max_plots elements
+        persistence = sorted(persistence, key=lambda life_time: life_time[1][1]-life_time[1][0], reverse=True)[:max_barcodes]
+
     (min_birth, max_death) = __min_birth_max_death(persistence)
     ind = 0
     delta = ((max_death - min_birth) / 10.0)
@@ -106,18 +131,40 @@ def plot_persistence_barcode(persistence, alpha=0.6):
     plt.title('Persistence barcode')
     # Ends plot on infinity value and starts a little bit before min_birth
     plt.axis([axis_start, infinity, 0, ind])
-    plt.show()
+    return plt
 
-def plot_persistence_diagram(persistence, alpha=0.6):
-    """This function plots the persistence diagram.
+def plot_persistence_diagram(persistence=[], persistence_file='', alpha=0.6, band_boot=0., max_plots=0):
+    """This function plots the persistence diagram with an optional confidence band.
 
     :param persistence: The persistence to plot.
     :type persistence: list of tuples(dimension, tuple(birth, death)).
+    :param persistence_file: A persistence file style name (reset persistence if both are set).
+    :type persistence_file: string
     :param alpha: alpha value in [0.0, 1.0] for points and horizontal infinity line (default is 0.6).
     :type alpha: float.
-    :returns: plot -- An diagram plot of persistence.
+    :param band_boot: bootstrap band (not displayed if :math:`\leq` 0.)
+    :type band_boot: float.
+    :param max_plots: number of maximal plots to be displayed
+    :type max_plots: int.
+    :returns: plot -- A diagram plot of persistence.
     """
-    (min_birth, max_death) = __min_birth_max_death(persistence)
+    if persistence_file is not '':
+        if os.path.isfile(persistence_file):
+            # Reset persistence
+            persistence = []
+            diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file)
+            for key in diag.keys():
+                for persistence_interval in diag[key]:
+                    persistence.append((key, persistence_interval))
+        else:
+            print("file " + persistence_file + " not found.")
+            return None
+
+    if max_plots > 0 and max_plots < len(persistence):
+        # Sort by life time, then takes only the max_plots elements
+        persistence = sorted(persistence, key=lambda life_time: life_time[1][1]-life_time[1][0], reverse=True)[:max_plots]
+
+    (min_birth, max_death) = __min_birth_max_death(persistence, band_boot)
     ind = 0
     delta = ((max_death - min_birth) / 10.0)
     # Replace infinity values with max_death + delta for diagram to be more
@@ -131,6 +178,9 @@ def plot_persistence_diagram(persistence, alpha=0.6):
     plt.plot(x, x, color='k', linewidth=1.0)
     plt.plot(x, [infinity] * len(x), linewidth=1.0, color='k', alpha=alpha)
     plt.text(axis_start, infinity, r'$\infty$', color='k', alpha=alpha)
+    # bootstrap band
+    if band_boot > 0.:
+        plt.fill_between(x, x, x+band_boot, alpha=alpha, facecolor='red')
 
     # Draw points in loop
     for interval in reversed(persistence):
@@ -149,4 +199,4 @@ def plot_persistence_diagram(persistence, alpha=0.6):
     plt.ylabel('Death')
     # Ends plot on infinity value and starts a little bit before min_birth
     plt.axis([axis_start, infinity, axis_start, infinity + delta])
-    plt.show()
+    return plt
diff --git a/src/cython/cython/reader_utils.pyx b/src/cython/cython/reader_utils.pyx
new file mode 100644
index 00000000..3a17c5a0
--- /dev/null
+++ b/src/cython/cython/reader_utils.pyx
@@ -0,0 +1,95 @@
+from cython cimport numeric
+from libcpp.vector cimport vector
+from libcpp.string cimport string
+from libcpp.map cimport map
+from libcpp.pair cimport pair
+import os
+
+"""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) 2017 INRIA
+
+   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/>.
+"""
+
+__author__ = "Vincent Rouvreau"
+__copyright__ = "Copyright (C) 2017 INRIA"
+__license__ = "GPL v3"
+
+cdef extern from "Reader_utils_interface.h" namespace "Gudhi":
+    vector[vector[double]] read_matrix_from_csv_file(string off_file, char separator)
+    map[int, vector[pair[double, double]]] read_pers_intervals_grouped_by_dimension(string filename)
+    vector[pair[double, double]] read_pers_intervals_in_dimension(string filename, int only_this_dim)
+
+def read_lower_triangular_matrix_from_csv_file(csv_file='', separator=';'):
+    """Read lower triangular matrix from a CSV style file.
+
+    :param csv_file: A CSV file style name.
+    :type csv_file: string
+    :param separator: The value separator in the CSV file. Default value is ';'
+    :type separator: char
+
+    :returns:  The lower triangular matrix.
+    :rtype: vector[vector[double]]
+    """
+    if csv_file is not '':
+        if os.path.isfile(csv_file):
+            return read_matrix_from_csv_file(str.encode(csv_file), ord(separator[0]))
+    print("file " + csv_file + " not set or not found.")
+    return []
+
+def read_persistence_intervals_grouped_by_dimension(persistence_file=''):
+    """Reads a file containing persistence intervals.
+    Each line might contain 2, 3 or 4 values: [[field] dimension] birth death
+    The return value is an `map[dim, vector[pair[birth, death]]]`
+    where `dim` is an `int`, `birth` a `double`, and `death` a `double`.
+    Note: the function does not check that birth <= death.
+
+    :param persistence_file: A persistence file style name.
+    :type persistence_file: string
+
+    :returns:  The persistence pairs grouped by dimension.
+    :rtype: map[int, vector[pair[double, double]]]
+    """
+    if persistence_file is not '':
+        if os.path.isfile(persistence_file):
+            return read_pers_intervals_grouped_by_dimension(str.encode(persistence_file))
+    print("file " + persistence_file + " not set or not found.")
+    return []
+
+def read_persistence_intervals_in_dimension(persistence_file='', only_this_dim=-1):
+    """Reads a file containing persistence intervals.
+    Each line might contain 2, 3 or 4 values: [[field] dimension] birth death
+    If `only_this_dim` = -1, dimension is ignored and all lines are returned.
+    If `only_this_dim` is >= 0, only the lines where dimension = `only_this_dim`
+    (or where dimension is not specified) are returned.
+    The return value is an `vector[pair[birth, death]]`
+    where `birth` a `double`, and `death` a `double`.
+    Note: the function does not check that birth <= death.
+
+    :param persistence_file: A persistence file style name.
+    :type persistence_file: string
+
+    :returns:  The persistence pairs grouped by dimension.
+    :rtype: map[int, vector[pair[double, double]]]
+    """
+    if persistence_file is not '':
+        if os.path.isfile(persistence_file):
+            return read_pers_intervals_in_dimension(str.encode(persistence_file), only_this_dim)
+    print("file " + persistence_file + " not set or not found.")
+    return []
diff --git a/src/cython/cython/simplex_tree.pyx b/src/cython/cython/simplex_tree.pyx
index 9d40a8b5..47aa5311 100644
--- a/src/cython/cython/simplex_tree.pyx
+++ b/src/cython/cython/simplex_tree.pyx
@@ -2,6 +2,7 @@ from cython cimport numeric
 from libcpp.vector cimport vector
 from libcpp.utility cimport pair
 from libcpp cimport bool
+from libcpp.string cimport string
 
 """This file is part of the Gudhi Library. The Gudhi library
    (Geometric Understanding in Higher Dimensions) is a generic C++
@@ -61,6 +62,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
         vector[int] betti_numbers()
         vector[int] persistent_betti_numbers(double from_value, double to_value)
         vector[pair[double,double]] intervals_in_dimension(int dimension)
+        void write_output_diagram(string diagram_file_name)
 
 # SimplexTree python interface
 cdef class SimplexTree:
@@ -183,10 +185,10 @@ cdef class SimplexTree:
         :returns:  true if the simplex was found, false otherwise.
         :rtype:  bool
         """
-        cdef vector[int] complex
+        cdef vector[int] csimplex
         for i in simplex:
-            complex.push_back(i)
-        return self.thisptr.find_simplex(complex)
+            csimplex.push_back(i)
+        return self.thisptr.find_simplex(csimplex)
 
     def insert(self, simplex, filtration=0.0):
         """This function inserts the given N-simplex and its subfaces with the
@@ -200,10 +202,10 @@ cdef class SimplexTree:
         :returns:  true if the simplex was found, false otherwise.
         :rtype:  bool
         """
-        cdef vector[int] complex
+        cdef vector[int] csimplex
         for i in simplex:
-            complex.push_back(i)
-        return self.thisptr.insert_simplex_and_subfaces(complex,
+            csimplex.push_back(i)
+        return self.thisptr.insert_simplex_and_subfaces(csimplex,
                                                         <double>filtration)
 
     def get_filtration(self):
@@ -232,35 +234,35 @@ cdef class SimplexTree:
         :returns:  The (simplices of the) skeleton of a maximum dimension.
         :rtype:  list of tuples(simplex, filtration)
         """
-        cdef vector[pair[vector[int], double]] skeletons \
+        cdef vector[pair[vector[int], double]] skeleton \
             = self.thisptr.get_skeleton(<int>dimension)
         ct = []
-        for filtered_complex in skeletons:
+        for filtered_simplex in skeleton:
             v = []
-            for vertex in filtered_complex.first:
+            for vertex in filtered_simplex.first:
                 v.append(vertex)
-            ct.append((v, filtered_complex.second))
+            ct.append((v, filtered_simplex.second))
         return ct
 
     def get_star(self, simplex):
-        """This function returns the stars of a given N-simplex.
+        """This function returns the star of a given N-simplex.
 
         :param simplex: The N-simplex, represented by a list of vertex.
         :type simplex: list of int.
         :returns:  The (simplices of the) star of a simplex.
         :rtype:  list of tuples(simplex, filtration)
         """
-        cdef vector[int] complex
+        cdef vector[int] csimplex
         for i in simplex:
-            complex.push_back(i)
-        cdef vector[pair[vector[int], double]] stars \
-            = self.thisptr.get_star(complex)
+            csimplex.push_back(i)
+        cdef vector[pair[vector[int], double]] star \
+            = self.thisptr.get_star(csimplex)
         ct = []
-        for filtered_complex in stars:
+        for filtered_simplex in star:
             v = []
-            for vertex in filtered_complex.first:
+            for vertex in filtered_simplex.first:
                 v.append(vertex)
-            ct.append((v, filtered_complex.second))
+            ct.append((v, filtered_simplex.second))
         return ct
 
     def get_cofaces(self, simplex, codimension):
@@ -275,17 +277,17 @@ cdef class SimplexTree:
         :returns:  The (simplices of the) cofaces of a simplex
         :rtype:  list of tuples(simplex, filtration)
         """
-        cdef vector[int] complex
+        cdef vector[int] csimplex
         for i in simplex:
-            complex.push_back(i)
+            csimplex.push_back(i)
         cdef vector[pair[vector[int], double]] cofaces \
-            = self.thisptr.get_cofaces(complex, <int>codimension)
+            = self.thisptr.get_cofaces(csimplex, <int>codimension)
         ct = []
-        for filtered_complex in cofaces:
+        for filtered_simplex in cofaces:
             v = []
-            for vertex in filtered_complex.first:
+            for vertex in filtered_simplex.first:
                 v.append(vertex)
-            ct.append((v, filtered_complex.second))
+            ct.append((v, filtered_simplex.second))
         return ct
 
     def remove_maximal_simplex(self, simplex):
@@ -385,7 +387,7 @@ cdef class SimplexTree:
         complex in a specific dimension.
 
         :param dimension: The specific dimension.
-        :type from_value: int.
+        :type dimension: int.
         :returns: The persistence intervals.
         :rtype:  list of pair of float
 
@@ -399,3 +401,22 @@ cdef class SimplexTree:
             print("intervals_in_dim function requires persistence function"
                   " to be launched first.")
         return intervals_result
+
+    def write_persistence_diagram(self, persistence_file=''):
+        """This function writes the persistence intervals of the simplicial
+        complex in a user given file name.
+
+        :param persistence_file: The specific dimension.
+        :type persistence_file: string.
+
+        :note: intervals_in_dim function requires persistence function to be
+            launched first.
+        """
+        if self.pcohptr != NULL:
+            if persistence_file != '':
+                self.pcohptr.write_output_diagram(str.encode(persistence_file))
+            else:
+                print("persistence_file must be specified")
+        else:
+            print("intervals_in_dim function requires persistence function"
+                  " to be launched first.")
diff --git a/src/cython/doc/_templates/layout.html b/src/cython/doc/_templates/layout.html
index b11c1236..243f33c6 100644
--- a/src/cython/doc/_templates/layout.html
+++ b/src/cython/doc/_templates/layout.html
@@ -65,6 +65,7 @@
             {#- old style sidebars: using blocks -- should be deprecated #}
             {%- block sidebartoc %}
 <h2><a href="index.html">GUDHI</a></h2>
+<h2><a href="fileformats.html">File formats</a></h2>
 <h2><a href="installation.html">GUDHI installation</a></h2>
 <h2><a href="citation.html">Acknowledging the GUDHI library</a></h2>
 <h2><a href="genindex.html">Index</a></h2>
diff --git a/src/cython/doc/alpha_complex_user.rst b/src/cython/doc/alpha_complex_user.rst
index e8268ef1..ff62390b 100644
--- a/src/cython/doc/alpha_complex_user.rst
+++ b/src/cython/doc/alpha_complex_user.rst
@@ -142,7 +142,7 @@ Prune above given filtration value
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 The simplex tree is pruned from the given maximum alpha squared value (cf. `Simplex_tree::prune_above_filtration()`
-int he `C++ version <http://gudhi.gforge.inria.fr/doc/latest/index.html>`_).
+in the `C++ version <http://gudhi.gforge.inria.fr/doc/latest/index.html>`_).
 In the following example, the value is given by the user as argument of the program.
 
 
diff --git a/src/cython/doc/conf.py b/src/cython/doc/conf.py
index 42bfd59c..072f7a94 100755
--- a/src/cython/doc/conf.py
+++ b/src/cython/doc/conf.py
@@ -69,7 +69,7 @@ copyright = u'2016, GUDHI Editorial Board'
 # The short X.Y version.
 version = '2.0'
 # The full version, including alpha/beta/rc tags.
-release = '2.0.0'
+release = '2.0.1-rc1'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/src/cython/doc/fileformats.rst b/src/cython/doc/fileformats.rst
new file mode 100644
index 00000000..156ef4e4
--- /dev/null
+++ b/src/cython/doc/fileformats.rst
@@ -0,0 +1,33 @@
+File formats
+############
+
+Persistence Diagram
+*******************
+
+Such a file, whose extension is usually ``.pers``, contains a list of
+persistence intervals.
+
+Lines starting with ``#`` are ignored (comments).
+
+Other lines might contain 2, 3 or 4 values (the number of values on each line
+must be the same for all lines)::
+
+    [[field] dimension] birth death
+
+Here is a simple sample file::
+
+    # Persistence diagram example
+    2 2.7 3.7
+    2 9.6 14.
+    # Some comments
+    3 34.2 34.974
+    4 3. inf
+
+Other sample files can be found in the data/persistence_diagram folder.
+
+Such files can be generated with
+:meth:`gudhi.SimplexTree.write_persistence_diagram`, read with
+:meth:`gudhi.read_persistence_intervals_grouped_by_dimension`, or
+:meth:`gudhi.read_persistence_intervals_in_dimension` and displayed with
+:meth:`gudhi.plot_persistence_barcode` or
+:meth:`gudhi.plot_persistence_diagram`.
diff --git a/src/cython/doc/persistence_graphical_tools_user.rst b/src/cython/doc/persistence_graphical_tools_user.rst
index cae18323..13198162 100644
--- a/src/cython/doc/persistence_graphical_tools_user.rst
+++ b/src/cython/doc/persistence_graphical_tools_user.rst
@@ -14,12 +14,14 @@ This function is useful to show the color palette values of dimension:
 .. testcode::
 
     import gudhi
-    gudhi.show_palette_values(alpha=1.0)
+    plt = gudhi.show_palette_values(alpha=1.0)
+    plt.show()
 
 .. plot::
 
     import gudhi
-    gudhi.show_palette_values(alpha=1.0)
+    plt = gudhi.show_palette_values(alpha=1.0)
+    plt.show()
 
 Show persistence as a barcode
 -----------------------------
@@ -32,7 +34,8 @@ This function can display the persistence result as a barcode:
     
     periodic_cc = gudhi.PeriodicCubicalComplex(perseus_file='3d_torus.txt')
     diag = periodic_cc.persistence()
-    gudhi.plot_persistence_barcode(diag)
+    plt = gudhi.plot_persistence_barcode(diag)
+    plt.show()
 
 .. plot::
 
@@ -40,7 +43,8 @@ This function can display the persistence result as a barcode:
 
     periodic_cc = gudhi.PeriodicCubicalComplex(perseus_file='3d_torus.txt')
     diag = periodic_cc.persistence()
-    gudhi.plot_persistence_barcode(diag)
+    plt = gudhi.plot_persistence_barcode(diag)
+    plt.show()
 
 Show persistence as a diagram
 -----------------------------
@@ -51,16 +55,18 @@ This function can display the persistence result as a diagram:
 
     import gudhi
     
-    rips_complex = gudhi.RipsComplex(off_file='tore3D_300.off', max_edge_length=2.0)
+    rips_complex = gudhi.RipsComplex(off_file='tore3D_1307.off', max_edge_length=0.2)
     simplex_tree = rips_complex.create_simplex_tree(max_dimension=3)
     diag = simplex_tree.persistence()
-    gudhi.plot_persistence_diagram(diag)
+    plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13)
+    plt.show()
 
 .. plot::
 
     import gudhi
 
-    rips_complex = gudhi.RipsComplex(off_file='tore3D_300.off', max_edge_length=2.0)
+    rips_complex = gudhi.RipsComplex(off_file='tore3D_1307.off', max_edge_length=0.2)
     simplex_tree = rips_complex.create_simplex_tree(max_dimension=3)
     diag = simplex_tree.persistence()
-    gudhi.plot_persistence_diagram(diag)
+    plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13)
+    plt.show()
diff --git a/src/cython/doc/pyplots/barcode_persistence.py b/src/cython/doc/pyplots/barcode_persistence.py
index c06ac5a7..9cd3149d 100755
--- a/src/cython/doc/pyplots/barcode_persistence.py
+++ b/src/cython/doc/pyplots/barcode_persistence.py
@@ -1,5 +1,6 @@
 import gudhi
 
-periodic_cc = gudhi.PeriodicCubicalComplex(perseus_file='../3d_torus.txt')
+periodic_cc = gudhi.PeriodicCubicalComplex(perseus_file='3d_torus.txt')
 diag = periodic_cc.persistence()
-gudhi.plot_persistence_barcode(diag)
+plt = gudhi.plot_persistence_barcode(diag)
+plt.show()
diff --git a/src/cython/doc/pyplots/diagram_persistence.py b/src/cython/doc/pyplots/diagram_persistence.py
index b4714fe3..30661965 100755
--- a/src/cython/doc/pyplots/diagram_persistence.py
+++ b/src/cython/doc/pyplots/diagram_persistence.py
@@ -1,5 +1,7 @@
 import gudhi
 
-alpha_complex = gudhi.AlphaComplex(off_file='../tore3D_300.off')
-diag = alpha_complex.persistence()
-gudhi.plot_persistence_diagram(diag)
+rips_complex = gudhi.RipsComplex(off_file='tore3D_1307.off', max_edge_length=0.2)
+simplex_tree = rips_complex.create_simplex_tree(max_dimension=3)
+diag = simplex_tree.persistence()
+plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13)
+plt.show()
diff --git a/src/cython/doc/pyplots/show_palette_values.py b/src/cython/doc/pyplots/show_palette_values.py
index e72a55fd..fdf9645f 100755
--- a/src/cython/doc/pyplots/show_palette_values.py
+++ b/src/cython/doc/pyplots/show_palette_values.py
@@ -1,2 +1,3 @@
 import gudhi
-gudhi.show_palette_values(alpha=1.0)
+plt = gudhi.show_palette_values(alpha=1.0)
+plt.show()
diff --git a/src/cython/doc/python3-sphinx-build b/src/cython/doc/python3-sphinx-build.in
index 44b94169..c97965f5 100755
--- a/src/cython/doc/python3-sphinx-build
+++ b/src/cython/doc/python3-sphinx-build.in
@@ -1,4 +1,4 @@
-#!/usr/bin/python3
+#!@PYTHON_EXECUTABLE@
 
 """
 Emulate sphinx-build for python3
diff --git a/src/cython/doc/reader_utils_ref.rst b/src/cython/doc/reader_utils_ref.rst
new file mode 100644
index 00000000..9c1ea6fc
--- /dev/null
+++ b/src/cython/doc/reader_utils_ref.rst
@@ -0,0 +1,11 @@
+=============================
+Reader utils reference manual
+=============================
+
+.. autofunction:: gudhi.read_off
+
+.. autofunction:: gudhi.read_lower_triangular_matrix_from_csv_file
+
+.. autofunction:: gudhi.read_persistence_intervals_grouped_by_dimension
+
+.. autofunction:: gudhi.read_persistence_intervals_in_dimension
diff --git a/src/cython/example/alpha_complex_diagram_persistence_from_off_file_example.py b/src/cython/example/alpha_complex_diagram_persistence_from_off_file_example.py
index adedc7d2..b4487be4 100755
--- a/src/cython/example/alpha_complex_diagram_persistence_from_off_file_example.py
+++ b/src/cython/example/alpha_complex_diagram_persistence_from_off_file_example.py
@@ -38,6 +38,7 @@ parser = argparse.ArgumentParser(description='AlphaComplex creation from '
                                  'points from the given OFF file.')
 parser.add_argument("-f", "--file", type=str, required=True)
 parser.add_argument("-a", "--max_alpha_square", type=float, default=0.5)
+parser.add_argument("-b", "--band_boot", type=float, default=0.)
 parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams')
 
 args = parser.parse_args()
@@ -63,7 +64,8 @@ with open(args.file, 'r') as f:
         print(simplex_tree.betti_numbers())
     
         if args.no_diagram == False:
-            gudhi.plot_persistence_diagram(diag)
+            pplot = gudhi.plot_persistence_diagram(diag, band_boot=args.band_boot)
+            pplot.show()
     else:
         print(args.file, "is not a valid OFF file")
 
diff --git a/src/cython/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py b/src/cython/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py
index 2371c36c..e3f362dc 100755
--- a/src/cython/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py
+++ b/src/cython/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py
@@ -40,6 +40,7 @@ parser.add_argument("-f", "--file", type=str, required=True)
 parser.add_argument("-a", "--max_alpha_square", type=float, required=True)
 parser.add_argument("-n", "--number_of_landmarks", type=int, required=True)
 parser.add_argument("-d", "--limit_dimension", type=int, required=True)
+parser.add_argument("-b", "--band_boot", type=float, default=0.)
 parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams')
 
 args = parser.parse_args()
@@ -70,8 +71,8 @@ with open(args.file, 'r') as f:
         print(simplex_tree.betti_numbers())
 
         if args.no_diagram == False:
-            gudhi.plot_persistence_diagram(diag)
-
+            pplot = gudhi.plot_persistence_diagram(diag, band_boot=args.band_boot)
+            pplot.show()
     else:
         print(args.file, "is not a valid OFF file")
 
diff --git a/src/cython/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py b/src/cython/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py
index 5748aa8a..c236d992 100755
--- a/src/cython/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py
+++ b/src/cython/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py
@@ -40,6 +40,7 @@ parser.add_argument("-f", "--file", type=str, required=True)
 parser.add_argument("-a", "--max_alpha_square", type=float, required=True)
 parser.add_argument("-n", "--number_of_landmarks", type=int, required=True)
 parser.add_argument("-d", "--limit_dimension", type=int, required=True)
+parser.add_argument("-b", "--band_boot", type=float, default=0.)
 parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams')
 
 args = parser.parse_args()
@@ -70,8 +71,8 @@ with open(args.file, 'r') as f:
         print(simplex_tree.betti_numbers())
 
         if args.no_diagram == False:
-            gudhi.plot_persistence_diagram(diag)
-
+            pplot = gudhi.plot_persistence_diagram(diag, band_boot=args.band_boot)
+            pplot.show()
     else:
         print(args.file, "is not a valid OFF file")
 
diff --git a/src/cython/example/gudhi_graphical_tools_example.py b/src/cython/example/gudhi_graphical_tools_example.py
index bc3b16ec..ed87806b 100755
--- a/src/cython/example/gudhi_graphical_tools_example.py
+++ b/src/cython/example/gudhi_graphical_tools_example.py
@@ -44,4 +44,11 @@ gudhi.plot_persistence_barcode(persistence)
 print("#####################################################################")
 print("Show diagram persistence example")
 
-gudhi.plot_persistence_diagram(persistence)
+pplot = gudhi.plot_persistence_diagram(persistence)
+pplot.show()
+
+print("#####################################################################")
+print("Show diagram persistence example with a confidence band")
+
+pplot = gudhi.plot_persistence_diagram(persistence, band_boot=0.2)
+pplot.show()
diff --git a/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py
index 984dbf1b..3baebd17 100755
--- a/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py
+++ b/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py
@@ -39,6 +39,7 @@ parser = argparse.ArgumentParser(description='RipsComplex creation from '
 parser.add_argument("-f", "--file", type=str, required=True)
 parser.add_argument("-e", "--max_edge_length", type=float, default=0.5)
 parser.add_argument("-d", "--max_dimension", type=int, default=1)
+parser.add_argument("-b", "--band_boot", type=float, default=0.)
 parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams')
 
 args = parser.parse_args()
@@ -61,4 +62,5 @@ print("betti_numbers()=")
 print(simplex_tree.betti_numbers())
 
 if args.no_diagram == False:
-    gudhi.plot_persistence_diagram(diag)
+    pplot = gudhi.plot_persistence_diagram(diag, band_boot=args.band_boot)
+    pplot.show()
diff --git a/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py
index 4c21b98e..5951eedf 100755
--- a/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py
+++ b/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py
@@ -39,6 +39,7 @@ parser = argparse.ArgumentParser(description='RipsComplex creation from '
 parser.add_argument("-f", "--file", type=str, required=True)
 parser.add_argument("-e", "--max_edge_length", type=float, default=0.5)
 parser.add_argument("-d", "--max_dimension", type=int, default=1)
+parser.add_argument("-b", "--band_boot", type=float, default=0.)
 parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams')
 
 args = parser.parse_args()
@@ -64,7 +65,8 @@ with open(args.file, 'r') as f:
         print(simplex_tree.betti_numbers())
     
         if args.no_diagram == False:
-            gudhi.plot_persistence_diagram(diag)
+            pplot = gudhi.plot_persistence_diagram(diag, band_boot=args.band_boot)
+            pplot.show()
     else:
         print(args.file, "is not a valid OFF file")
 
diff --git a/src/cython/example/rips_persistence_diagram.py b/src/cython/example/rips_persistence_diagram.py
index 4e5cd2c8..9bfea41c 100755
--- a/src/cython/example/rips_persistence_diagram.py
+++ b/src/cython/example/rips_persistence_diagram.py
@@ -39,4 +39,5 @@ simplex_tree = rips.create_simplex_tree(max_dimension=1)
 diag = simplex_tree.persistence(homology_coeff_field=2, min_persistence=0)
 print("diag=", diag)
 
-gudhi.plot_persistence_diagram(diag)
+pplot = gudhi.plot_persistence_diagram(diag)
+pplot.show()
diff --git a/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py b/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py
index 4845eb47..6145e7f2 100755
--- a/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py
+++ b/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py
@@ -37,6 +37,7 @@ parser = argparse.ArgumentParser(description='TangentialComplex creation from '
                                  '- Constructs a tangential complex with the '
                                  'points from the given OFF file')
 parser.add_argument("-f", "--file", type=str, required=True)
+parser.add_argument("-b", "--band_boot", type=float, default=0.)
 parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams')
 
 args = parser.parse_args()
@@ -59,7 +60,8 @@ with open(args.file, 'r') as f:
         print(st.betti_numbers())
     
         if args.no_diagram == False:
-            gudhi.plot_persistence_diagram(diag)
+            pplot = gudhi.plot_persistence_diagram(diag, band_boot=args.band_boot)
+            pplot.show()
     else:
         print(args.file, "is not a valid OFF file")
 
diff --git a/src/cython/gudhi.pyx.in b/src/cython/gudhi.pyx.in
index 34d7c3b5..ed2d28cc 100644
--- a/src/cython/gudhi.pyx.in
+++ b/src/cython/gudhi.pyx.in
@@ -30,6 +30,7 @@ include "cython/rips_complex.pyx"
 include "cython/cubical_complex.pyx"
 include "cython/periodic_cubical_complex.pyx"
 include "cython/persistence_graphical_tools.py"
+include "cython/reader_utils.pyx"
 include "cython/witness_complex.pyx"
 include "cython/strong_witness_complex.pyx"
 @GUDHI_CYTHON_ALPHA_COMPLEX@
diff --git a/src/cython/include/Reader_utils_interface.h b/src/cython/include/Reader_utils_interface.h
new file mode 100644
index 00000000..b87b6cca
--- /dev/null
+++ b/src/cython/include/Reader_utils_interface.h
@@ -0,0 +1,54 @@
+/*    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) 2017 INRIA
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef INCLUDE_READER_UTILS_INTERFACE_H_
+#define INCLUDE_READER_UTILS_INTERFACE_H_
+
+#include <gudhi/reader_utils.h>
+
+#include <iostream>
+#include <vector>
+#include <string>
+
+namespace Gudhi {
+
+// Redefine functions with a different name in order the original name can be used in the Python version.
+std::vector<std::vector<double>> read_matrix_from_csv_file(const std::string& filename,
+                                                                const char separator = ';') {
+  return read_lower_triangular_matrix_from_csv_file<double>(filename, separator);
+}
+
+inline std::map<int, std::vector<std::pair<double, double>>>
+    read_pers_intervals_grouped_by_dimension(std::string const& filename) {
+  return read_persistence_intervals_grouped_by_dimension(filename);
+}
+
+inline std::vector<std::pair<double, double>>
+    read_pers_intervals_in_dimension(std::string const& filename, int only_this_dim = -1) {
+  return read_persistence_intervals_in_dimension(filename, only_this_dim);
+}
+
+
+}  // namespace Gudhi
+
+
+#endif  // INCLUDE_READER_UTILS_INTERFACE_H_
diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h
index 6d813f4a..d06ee4bd 100644
--- a/src/cython/include/Rips_complex_interface.h
+++ b/src/cython/include/Rips_complex_interface.h
@@ -66,7 +66,7 @@ class Rips_complex_interface {
     } else {
       // Rips construction where values is a distance matrix
       Distance_matrix distances =
-          read_lower_triangular_matrix_from_csv_file<Simplex_tree_interface<>::Filtration_value>(file_name);
+          Gudhi::read_lower_triangular_matrix_from_csv_file<Simplex_tree_interface<>::Filtration_value>(file_name);
       rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(distances, threshold);
     }
   }
diff --git a/src/cython/include/Simplex_tree_interface.h b/src/cython/include/Simplex_tree_interface.h
index 45ce1916..09e7e992 100644
--- a/src/cython/include/Simplex_tree_interface.h
+++ b/src/cython/include/Simplex_tree_interface.h
@@ -70,14 +70,14 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
   }
 
   // Do not interface this function, only used in strong witness interface for complex creation
-  bool insert_simplex(const std::vector<std::size_t>& complex, Filtration_value filtration = 0) {
-    Insertion_result result = Base::insert_simplex(complex, filtration);
+  bool insert_simplex(const std::vector<std::size_t>& simplex, Filtration_value filtration = 0) {
+    Insertion_result result = Base::insert_simplex(simplex, filtration);
     return (result.second);
   }
 
   // Do not interface this function, only used in strong witness interface for complex creation
-  bool insert_simplex_and_subfaces(const std::vector<std::size_t>& complex, Filtration_value filtration = 0) {
-    Insertion_result result = Base::insert_simplex_and_subfaces(complex, filtration);
+  bool insert_simplex_and_subfaces(const std::vector<std::size_t>& simplex, Filtration_value filtration = 0) {
+    Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration);
     return (result.second);
   }
 
diff --git a/src/cython/test/test_cubical_complex.py b/src/cython/test/test_cubical_complex.py
index 2e281ee4..9a365823 100755
--- a/src/cython/test/test_cubical_complex.py
+++ b/src/cython/test/test_cubical_complex.py
@@ -67,7 +67,7 @@ def test_dimension_constructor():
                          top_dimensional_cells = [1,2,3,4,5,6,7,8,9])
     assert cub.__is_defined() == True
     assert cub.__is_persistence_defined() == False
-    assert cub.persistence() == [(1, (0.0, 100.0)), (0, (0.0, 1.8446744073709552e+19))]
+    assert cub.persistence() == [(1, (0.0, 100.0)), (0, (0.0, float('inf')))]
     assert cub.__is_persistence_defined() == True
     assert cub.betti_numbers() == [1, 0]
     assert cub.persistent_betti_numbers(0, 1000) == [0, 0]
@@ -80,7 +80,7 @@ def test_dimension_constructor():
     cub = CubicalComplex(perseus_file='CubicalOneSphere.txt')
     assert cub.__is_defined() == True
     assert cub.__is_persistence_defined() == False
-    assert cub.persistence() == [(1, (0.0, 100.0)), (0, (0.0, 1.8446744073709552e+19))]
+    assert cub.persistence() == [(1, (0.0, 100.0)), (0, (0.0, float('inf')))]
     assert cub.__is_persistence_defined() == True
     assert cub.betti_numbers() == [1, 0, 0]
     assert cub.persistent_betti_numbers(0, 1000) == [1, 0, 0]
diff --git a/src/cython/test/test_reader_utils.py b/src/cython/test/test_reader_utils.py
new file mode 100755
index 00000000..25591fb3
--- /dev/null
+++ b/src/cython/test/test_reader_utils.py
@@ -0,0 +1,88 @@
+import gudhi
+
+"""This file is part of the Gudhi Library. The Gudhi library
+   (Geometric Understanding in Higher Dimensions) is a generic C++
+   library for computational topology.
+
+   Author(s):       Vincent Rouvreau
+
+   Copyright (C) 2017 INRIA
+
+   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/>.
+"""
+
+__author__ = "Vincent Rouvreau"
+__copyright__ = "Copyright (C) 2017 INRIA"
+__license__ = "GPL v3"
+
+
+def test_non_existing_csv_file():
+    # Try to open a non existing file
+    matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file='pouetpouettralala.toubiloubabdou')
+    assert matrix == []
+
+def test_full_square_distance_matrix_csv_file():
+    # Create test file
+    test_file = open('full_square_distance_matrix.csv', 'w')
+    test_file.write('0;1;2;3;\n1;0;4;5;\n2;4;0;6;\n3;5;6;0;')
+    test_file.close()
+    matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file="full_square_distance_matrix.csv")
+    assert matrix == [[], [1.0], [2.0, 4.0], [3.0, 5.0, 6.0]]
+
+def test_lower_triangular_distance_matrix_csv_file():
+    # Create test file
+    test_file = open('lower_triangular_distance_matrix.csv', 'w')
+    test_file.write('\n1,\n2,3,\n4,5,6,\n7,8,9,10,')
+    test_file.close()
+    matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file="lower_triangular_distance_matrix.csv", separator=",")
+    assert matrix == [[], [1.0], [2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0, 10.0]]
+
+def test_non_existing_persistence_file():
+    # Try to open a non existing file
+    persistence = gudhi.read_persistence_intervals_grouped_by_dimension(persistence_file='pouetpouettralala.toubiloubabdou')
+    assert persistence == []
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='pouetpouettralala.toubiloubabdou', only_this_dim=1)
+    assert persistence == []
+
+def test_read_persistence_intervals_without_dimension():
+    # Create test file
+    test_file = open('persistence_intervals_without_dimension.pers', 'w')
+    test_file.write('# Simple persistence diagram without dimension\n2.7 3.7\n9.6 14.\n34.2 34.974\n3. inf')
+    test_file.close()
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='persistence_intervals_without_dimension.pers')
+    assert persistence == [(2.7, 3.7), (9.6, 14.), (34.2, 34.974), (3., float('Inf'))]
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='persistence_intervals_without_dimension.pers', only_this_dim=0)
+    assert persistence == []
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='persistence_intervals_without_dimension.pers', only_this_dim=1)
+    assert persistence == []
+    persistence = gudhi.read_persistence_intervals_grouped_by_dimension(persistence_file='persistence_intervals_without_dimension.pers')
+    assert persistence == {-1: [(2.7, 3.7), (9.6, 14.0), (34.2, 34.974), (3.0, float('Inf'))]}
+
+def test_read_persistence_intervals_with_dimension():
+    # Create test file
+    test_file = open('persistence_intervals_with_dimension.pers', 'w')
+    test_file.write('# Simple persistence diagram with dimension\n0 2.7 3.7\n1 9.6 14.\n3 34.2 34.974\n1 3. inf')
+    test_file.close()
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='persistence_intervals_with_dimension.pers')
+    assert persistence == [(2.7, 3.7), (9.6, 14.), (34.2, 34.974), (3., float('Inf'))]
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='persistence_intervals_with_dimension.pers', only_this_dim=0)
+    assert persistence == [(2.7, 3.7)]
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='persistence_intervals_with_dimension.pers', only_this_dim=1)
+    assert persistence == [(9.6, 14.), (3., float('Inf'))]
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='persistence_intervals_with_dimension.pers', only_this_dim=2)
+    assert persistence == []
+    persistence = gudhi.read_persistence_intervals_in_dimension(persistence_file='persistence_intervals_with_dimension.pers', only_this_dim=3)
+    assert persistence == [(34.2, 34.974)]
+    persistence = gudhi.read_persistence_intervals_grouped_by_dimension(persistence_file='persistence_intervals_with_dimension.pers')
+    assert persistence == {0: [(2.7, 3.7)], 1: [(9.6, 14.0), (3.0, float('Inf'))], 3: [(34.2, 34.974)]}