From cea821f9ca34c270a5ccc047342c2c21ae79a6c0 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Fri, 28 Aug 2020 17:42:12 +0200
Subject: A prototype to fix #364

---
 src/python/gudhi/simplex_tree.pxd           |  3 ++
 src/python/gudhi/simplex_tree.pyx           | 21 ++++++++++++++
 src/python/include/Simplex_tree_interface.h | 12 ++++++++
 src/python/test/test_simplex_tree.py        | 44 +++++++++++++++++++++++++++++
 4 files changed, 80 insertions(+)

diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 75e94e0b..44533d7f 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -66,6 +66,9 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_end() nogil
         Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil
         Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil
+        #
+        ctypedef bool (*blocker_func)(vector[int], void *user_data)
+        void expansion_with_blockers_callback(int dimension, blocker_func user_func, void *user_data)
 
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>":
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index dfb1d985..e297cfbd 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -17,6 +17,9 @@ __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
 __license__ = "MIT"
 
+cdef bool callback(vector[int] simplex, void *blocker_func):
+    return (<object>blocker_func)(simplex)
+
 # SimplexTree python interface
 cdef class SimplexTree:
     """The simplex tree is an efficient and flexible data structure for
@@ -409,6 +412,24 @@ cdef class SimplexTree:
         persistence_result = self.pcohptr.get_persistence()
         return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result, min_persistence)
 
+    def expansion_with_blocker(self, max_dim, blocker_func):
+        """Expands the Simplex_tree containing only its one skeleton
+        until dimension max_dim.
+
+        The expanded simplicial complex until dimension :math:`d`
+        attached to a graph :math:`G` is the maximal simplicial complex of
+        dimension at most :math:`d` admitting the graph :math:`G` as
+        :math:`1`-skeleton.
+        The filtration value assigned to a simplex is the maximal filtration
+        value of one of its edges.
+
+        The Simplex_tree must contain no simplex of dimension bigger than
+        1 when calling the method.
+
+        :param max_dim: The maximal dimension.
+        :type max_dim: int
+        """
+        self.get_ptr().expansion_with_blockers_callback(max_dim, callback, <void*>blocker_func)
 
     def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False):
         """This function computes and returns the persistence of the simplicial complex.
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index e288a8cf..9f92b349 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -39,6 +39,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
   using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator;
   using Complex_simplex_iterator = typename Base::Complex_simplex_iterator;
   using Extended_filtration_data = typename Base::Extended_filtration_data;
+  typedef bool (*blocker_func)(Simplex simplex, void *user_data);
 
  public:
 
@@ -188,6 +189,17 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     return collapsed_stree_ptr;
   }
 
+  void expansion_with_blockers_callback(int dimension, blocker_func user_func, void *user_data) {
+    Base::expansion_with_blockers(dimension, [&](Simplex_handle sh){
+      Simplex simplex;
+      for (auto vertex : Base::simplex_vertex_range(sh)) {
+        simplex.insert(simplex.begin(), vertex);
+      }
+      return user_func(simplex, user_data);
+    });
+    Base::clear_filtration();
+  }
+
   // Iterator over the simplex tree
   Complex_simplex_iterator get_simplices_iterator_begin() {
     // this specific case works because the range is just a pair of iterators - won't work if range was a vector
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 83be0602..7aad8259 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -358,3 +358,47 @@ def test_collapse_edges():
     assert st.find([1, 3]) == False
     for simplex in st.get_skeleton(0): 
         assert simplex[1] == 1. 
+
+def blocker(simplex):
+    try:
+        # Block all simplices that countains vertex 6
+        simplex.index(6)
+        print(simplex, ' is blocked')
+        return True
+    except ValueError:
+        print(simplex, ' is accepted')
+        return False
+
+def test_expansion_with_blocker():
+    st=SimplexTree()
+    st.insert([0,1],0)
+    st.insert([0,2],1)
+    st.insert([0,3],2)
+    st.insert([1,2],3)
+    st.insert([1,3],4)
+    st.insert([2,3],5)
+    st.insert([2,4],6)
+    st.insert([3,6],7)
+    st.insert([4,5],8)
+    st.insert([4,6],9)
+    st.insert([5,6],10)
+    st.insert([6],10)
+
+    st.expansion_with_blocker(2, blocker)
+    assert st.num_simplices() == 22
+    assert st.dimension() == 2
+    assert st.find([4,5,6]) == False
+    assert st.filtration([0,1,2]) == 3.
+    assert st.filtration([0,1,3]) == 4.
+    assert st.filtration([0,2,3]) == 5.
+    assert st.filtration([1,2,3]) == 5.
+
+    st.expansion_with_blocker(3, blocker)
+    assert st.num_simplices() == 23
+    assert st.dimension() == 3
+    assert st.find([4,5,6]) == False
+    assert st.filtration([0,1,2]) == 3.
+    assert st.filtration([0,1,3]) == 4.
+    assert st.filtration([0,2,3]) == 5.
+    assert st.filtration([1,2,3]) == 5.
+    assert st.filtration([0,1,2,3]) == 5.
-- 
cgit v1.2.3


From c2eb0484191f89fcbe40bc4ab04943eb808f12a9 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 1 Sep 2020 14:51:25 +0200
Subject: expansion with blocker and how to modify filtration value

---
 src/python/gudhi/simplex_tree.pxd    |  2 +-
 src/python/gudhi/simplex_tree.pyx    | 23 +++++++++---------
 src/python/test/test_simplex_tree.py | 46 +++++++++++++++++++++---------------
 3 files changed, 40 insertions(+), 31 deletions(-)

diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 44533d7f..80c6ffca 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -66,7 +66,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_end() nogil
         Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil
         Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil
-        #
+        # Expansion with blockers
         ctypedef bool (*blocker_func)(vector[int], void *user_data)
         void expansion_with_blockers_callback(int dimension, blocker_func user_func, void *user_data)
 
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index e297cfbd..debe92c0 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -413,21 +413,22 @@ cdef class SimplexTree:
         return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result, min_persistence)
 
     def expansion_with_blocker(self, max_dim, blocker_func):
-        """Expands the Simplex_tree containing only its one skeleton
-        until dimension max_dim.
+        """Expands the Simplex_tree containing only a graph. Simplices corresponding to cliques in the graph are added
+        incrementally, faces before cofaces, unless the simplex has dimension larger than `max_dim` or `blocker_func`
+        returns `True` for this simplex.
 
-        The expanded simplicial complex until dimension :math:`d`
-        attached to a graph :math:`G` is the maximal simplicial complex of
-        dimension at most :math:`d` admitting the graph :math:`G` as
-        :math:`1`-skeleton.
-        The filtration value assigned to a simplex is the maximal filtration
-        value of one of its edges.
+        The function identifies a candidate simplex whose faces are all already in the complex, inserts it with a
+        filtration value corresponding to the maximum of the filtration values of the faces, then calls `blocker_func`
+        with this new simplex (represented as a list of int). If `blocker_func` returns `True`, the simplex is removed,
+        otherwise it is kept. The algorithm then proceeds with the next candidate.
 
-        The Simplex_tree must contain no simplex of dimension bigger than
-        1 when calling the method.
+        Note that you cannot update the filtration value of the simplex during the evaluation of `blocker_func`, as it
+        would segfault. 
 
-        :param max_dim: The maximal dimension.
+        :param max_dim: Expansion maximal dimension value.
         :type max_dim: int
+        :param blocker_func: Blocker oracle.
+        :type blocker_func: Its concept is `Boolean blocker_func(list of int)`
         """
         self.get_ptr().expansion_with_blockers_callback(max_dim, callback, <void*>blocker_func)
 
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 7aad8259..33b0ac99 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -359,16 +359,6 @@ def test_collapse_edges():
     for simplex in st.get_skeleton(0): 
         assert simplex[1] == 1. 
 
-def blocker(simplex):
-    try:
-        # Block all simplices that countains vertex 6
-        simplex.index(6)
-        print(simplex, ' is blocked')
-        return True
-    except ValueError:
-        print(simplex, ' is accepted')
-        return False
-
 def test_expansion_with_blocker():
     st=SimplexTree()
     st.insert([0,1],0)
@@ -384,21 +374,39 @@ def test_expansion_with_blocker():
     st.insert([5,6],10)
     st.insert([6],10)
 
+    # One cannot modify filtration inside blocker - segfault - Let's record accepted simplices
+    accepted_simp = []
+    def blocker(simplex):
+        try:
+            # Block all simplices that countains vertex 6
+            simplex.index(6)
+            print(simplex, ' is blocked')
+            return True
+        except ValueError:
+            print(simplex, ' is accepted')
+            accepted_simp.append(simplex)
+            return False
+
     st.expansion_with_blocker(2, blocker)
+    for simplex in accepted_simp:
+        st.assign_filtration(simplex, st.filtration(simplex) + 1.)
     assert st.num_simplices() == 22
     assert st.dimension() == 2
     assert st.find([4,5,6]) == False
-    assert st.filtration([0,1,2]) == 3.
-    assert st.filtration([0,1,3]) == 4.
-    assert st.filtration([0,2,3]) == 5.
-    assert st.filtration([1,2,3]) == 5.
+    assert st.filtration([0,1,2]) == 4.
+    assert st.filtration([0,1,3]) == 5.
+    assert st.filtration([0,2,3]) == 6.
+    assert st.filtration([1,2,3]) == 6.
 
+    accepted_simp = []
     st.expansion_with_blocker(3, blocker)
+    for simplex in accepted_simp:
+        st.assign_filtration(simplex, st.filtration(simplex) + 1.)
     assert st.num_simplices() == 23
     assert st.dimension() == 3
     assert st.find([4,5,6]) == False
-    assert st.filtration([0,1,2]) == 3.
-    assert st.filtration([0,1,3]) == 4.
-    assert st.filtration([0,2,3]) == 5.
-    assert st.filtration([1,2,3]) == 5.
-    assert st.filtration([0,1,2,3]) == 5.
+    assert st.filtration([0,1,2]) == 4.
+    assert st.filtration([0,1,3]) == 5.
+    assert st.filtration([0,2,3]) == 6.
+    assert st.filtration([1,2,3]) == 6.
+    assert st.filtration([0,1,2,3]) == 6.
-- 
cgit v1.2.3


From f4f5992e9025d055e35c86589923f1b0299f552e Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Thu, 22 Apr 2021 08:39:03 +0200
Subject: code review: rename block_func when it is a type. Type is a
 'callable'

---
 src/python/gudhi/simplex_tree.pxd           | 4 ++--
 src/python/gudhi/simplex_tree.pyx           | 2 +-
 src/python/include/Simplex_tree_interface.h | 4 ++--
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 80c6ffca..ff750af9 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -67,8 +67,8 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil
         Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil
         # Expansion with blockers
-        ctypedef bool (*blocker_func)(vector[int], void *user_data)
-        void expansion_with_blockers_callback(int dimension, blocker_func user_func, void *user_data)
+        ctypedef bool (*blocker_func_t)(vector[int], void *user_data)
+        void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data)
 
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>":
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index debe92c0..383d1949 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -428,7 +428,7 @@ cdef class SimplexTree:
         :param max_dim: Expansion maximal dimension value.
         :type max_dim: int
         :param blocker_func: Blocker oracle.
-        :type blocker_func: Its concept is `Boolean blocker_func(list of int)`
+        :type blocker_func: Callable[[List[int]], bool]
         """
         self.get_ptr().expansion_with_blockers_callback(max_dim, callback, <void*>blocker_func)
 
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index 9f92b349..a859bec0 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -39,7 +39,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
   using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator;
   using Complex_simplex_iterator = typename Base::Complex_simplex_iterator;
   using Extended_filtration_data = typename Base::Extended_filtration_data;
-  typedef bool (*blocker_func)(Simplex simplex, void *user_data);
+  typedef bool (*blocker_func_t)(Simplex simplex, void *user_data);
 
  public:
 
@@ -189,7 +189,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     return collapsed_stree_ptr;
   }
 
-  void expansion_with_blockers_callback(int dimension, blocker_func user_func, void *user_data) {
+  void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) {
     Base::expansion_with_blockers(dimension, [&](Simplex_handle sh){
       Simplex simplex;
       for (auto vertex : Base::simplex_vertex_range(sh)) {
-- 
cgit v1.2.3


From 16c421a072d2f094b0541e7ef95a50427d39fec3 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 27 Apr 2021 09:03:14 +0200
Subject: This test was forgotten

---
 src/Simplex_tree/test/CMakeLists.txt | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/Simplex_tree/test/CMakeLists.txt b/src/Simplex_tree/test/CMakeLists.txt
index cf2b0153..25b562e0 100644
--- a/src/Simplex_tree/test/CMakeLists.txt
+++ b/src/Simplex_tree/test/CMakeLists.txt
@@ -34,3 +34,9 @@ if (TBB_FOUND)
   target_link_libraries(Simplex_tree_make_filtration_non_decreasing_test_unit ${TBB_LIBRARIES})
 endif()
 gudhi_add_boost_test(Simplex_tree_make_filtration_non_decreasing_test_unit)
+
+add_executable ( Simplex_tree_graph_expansion_test_unit simplex_tree_graph_expansion_unit_test.cpp )
+if (TBB_FOUND)
+  target_link_libraries(Simplex_tree_graph_expansion_test_unit ${TBB_LIBRARIES})
+endif()
+gudhi_add_boost_test(Simplex_tree_graph_expansion_test_unit)
-- 
cgit v1.2.3


From 69168e8ed24165ab89ea1c57bc21dd994c93dd8e Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 27 Apr 2021 09:03:47 +0200
Subject: vector initialization improvement

---
 src/python/include/Simplex_tree_interface.h | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index a859bec0..dd8d715c 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -191,13 +191,9 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
 
   void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) {
     Base::expansion_with_blockers(dimension, [&](Simplex_handle sh){
-      Simplex simplex;
-      for (auto vertex : Base::simplex_vertex_range(sh)) {
-        simplex.insert(simplex.begin(), vertex);
-      }
+      Simplex simplex(Base::simplex_vertex_range(sh).begin(), Base::simplex_vertex_range(sh).end());
       return user_func(simplex, user_data);
     });
-    Base::clear_filtration();
   }
 
   // Iterator over the simplex tree
-- 
cgit v1.2.3


From 17da2842a3af3db3abec7c2ce03f77dad47ef5dc Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 15 Jun 2021 18:01:20 +0200
Subject: Remove deprecated max_barcodes in plot_persistence_barcode and
 max_plots in plot_persistence_diagram. Fix #453 and factorize code

---
 src/python/gudhi/persistence_graphical_tools.py | 63 ++++++++++++-------------
 1 file changed, 29 insertions(+), 34 deletions(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 848dc03e..460e2558 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -12,6 +12,7 @@ from os import path
 from math import isfinite
 import numpy as np
 from functools import lru_cache
+import warnings
 
 from gudhi.reader_utils import read_persistence_intervals_in_dimension
 from gudhi.reader_utils import read_persistence_intervals_grouped_by_dimension
@@ -58,6 +59,26 @@ def _array_handler(a):
     else:
         return a
 
+def _limit_to_max_intervals(persistence, max_intervals, key):
+    '''This function returns truncated persistence if length is bigger than max_intervals.
+    :param persistence: Persistence intervals values list. Can be grouped by dimension or not.
+    :type persistence: an array of (dimension, array of (birth, death)) or an array of (birth, death).
+    :param max_intervals: maximal number of intervals to display.
+        Selected intervals are those with the longest life time. Set it
+        to 0 to see all. Default value is 1000.
+    :type max_intervals: int.
+    :param key: key function for sort algorithm.
+    :type key: function or lambda.
+    '''
+    if max_intervals > 0 and max_intervals < len(persistence):
+        warnings.warn('There are %s intervals given as input, whereas max_intervals is set to %s.' %
+                        (len(persistence), max_intervals)
+                     )
+        # Sort by life time, then takes only the max_intervals elements
+        return sorted(persistence, key = key, reverse = True)[:max_intervals]
+    else:
+        return persistence
+
 @lru_cache(maxsize=1)
 def _matplotlib_can_use_tex():
     """This function returns True if matplotlib can deal with LaTeX, False otherwise.
@@ -75,7 +96,6 @@ def plot_persistence_barcode(
     persistence_file="",
     alpha=0.6,
     max_intervals=1000,
-    max_barcodes=1000,
     inf_delta=0.1,
     legend=False,
     colormap=None,
@@ -142,18 +162,9 @@ def plot_persistence_barcode(
 
         persistence = _array_handler(persistence)
 
-        if max_barcodes != 1000:
-            print("Deprecated parameter. It has been replaced by max_intervals")
-            max_intervals = max_barcodes
-
-        if max_intervals > 0 and max_intervals < len(persistence):
-            # Sort by life time, then takes only the max_intervals elements
-            persistence = sorted(
-                persistence,
-                key=lambda life_time: life_time[1][1] - life_time[1][0],
-                reverse=True,
-            )[:max_intervals]
-        
+        persistence = _limit_to_max_intervals(persistence, max_intervals,
+                                              key = lambda life_time: life_time[1][1] - life_time[1][0])
+
         if colormap == None:
             colormap = plt.cm.Set1.colors
         if axes == None:
@@ -220,7 +231,6 @@ def plot_persistence_diagram(
     alpha=0.6,
     band=0.0,
     max_intervals=1000,
-    max_plots=1000,
     inf_delta=0.1,
     legend=False,
     colormap=None,
@@ -291,17 +301,8 @@ def plot_persistence_diagram(
 
         persistence = _array_handler(persistence)
 
-        if max_plots != 1000:
-            print("Deprecated parameter. It has been replaced by max_intervals")
-            max_intervals = max_plots
-
-        if max_intervals > 0 and max_intervals < len(persistence):
-            # Sort by life time, then takes only the max_intervals elements
-            persistence = sorted(
-                persistence,
-                key=lambda life_time: life_time[1][1] - life_time[1][0],
-                reverse=True,
-            )[:max_intervals]
+        persistence = _limit_to_max_intervals(persistence, max_intervals,
+                                              key = lambda life_time: life_time[1][1] - life_time[1][0])
 
         if colormap == None:
             colormap = plt.cm.Set1.colors
@@ -474,15 +475,9 @@ def plot_persistence_density(
             )
 
         persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])]
-        if max_intervals > 0 and max_intervals < len(persistence_dim):
-            # Sort by life time, then takes only the max_intervals elements
-            persistence_dim = np.array(
-                sorted(
-                    persistence_dim,
-                    key=lambda life_time: life_time[1] - life_time[0],
-                    reverse=True,
-                )[:max_intervals]
-            )
+
+        persistence_dim = np.array(_limit_to_max_intervals(persistence_dim, max_intervals,
+                                                           key = lambda life_time: life_time[1] - life_time[0]))
 
         # Set as numpy array birth and death (remove undefined values - inf and NaN)
         birth = persistence_dim[:, 0]
-- 
cgit v1.2.3


From 8587b70b268aaf7be2ab5b0e163d43f5587b0aac Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 15 Jun 2021 18:03:45 +0200
Subject: Fix MatplotlibDeprecationWarning: shading='flat' when X and Y have
 the same dimensions as C as it is deprecated since 3.3

---
 src/python/gudhi/persistence_graphical_tools.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 460e2558..460029fb 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -502,7 +502,7 @@ def plot_persistence_density(
         zi = k(np.vstack([xi.flatten(), yi.flatten()]))
 
         # Make the plot
-        img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap)
+        img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading='auto')
 
         if greyblock:
             axes.add_patch(mpatches.Polygon([[birth.min(), birth.min()], [death.max(), birth.min()], [death.max(), death.max()]], fill=True, color='lightgrey'))
-- 
cgit v1.2.3


From 2eb2f96248f09e5c793760cfb44ce7b1a40c1af4 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 15 Jun 2021 18:32:33 +0200
Subject: Ignore figure as not used on default subplot

---
 src/python/gudhi/persistence_graphical_tools.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 460029fb..b129b375 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -168,7 +168,7 @@ def plot_persistence_barcode(
         if colormap == None:
             colormap = plt.cm.Set1.colors
         if axes == None:
-            fig, axes = plt.subplots(1, 1)
+            _, axes = plt.subplots(1, 1)
 
         persistence = sorted(persistence, key=lambda birth: birth[1][0])
 
@@ -307,7 +307,7 @@ def plot_persistence_diagram(
         if colormap == None:
             colormap = plt.cm.Set1.colors
         if axes == None:
-            fig, axes = plt.subplots(1, 1)
+            _, axes = plt.subplots(1, 1)
 
         (min_birth, max_death) = __min_birth_max_death(persistence, band)
         delta = (max_death - min_birth) * inf_delta
@@ -487,7 +487,7 @@ def plot_persistence_density(
         if cmap is None:
             cmap = plt.cm.hot_r
         if axes == None:
-            fig, axes = plt.subplots(1, 1)
+            _, axes = plt.subplots(1, 1)
 
         # line display of equation : birth = death
         x = np.linspace(death.min(), birth.max(), 1000)
-- 
cgit v1.2.3


From 486b281c726cbb6110cfe3c63b3f225690bcd348 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Thu, 17 Jun 2021 11:01:18 +0200
Subject: Use float as it can be ambiguous. Do not work if we remove np.inf
 value

---
 src/python/doc/persistence_graphical_tools_user.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/doc/persistence_graphical_tools_user.rst b/src/python/doc/persistence_graphical_tools_user.rst
index d95b9d2b..e1d28c71 100644
--- a/src/python/doc/persistence_graphical_tools_user.rst
+++ b/src/python/doc/persistence_graphical_tools_user.rst
@@ -60,7 +60,7 @@ of shape (N x 2) encoding a persistence diagram (in a given dimension).
     import matplotlib.pyplot as plt
     import gudhi
     import numpy as np
-    d = np.array([[0, 1], [1, 2], [1, np.inf]])
+    d = np.array([[0., 1.], [1., 2.], [1., np.inf]])
     gudhi.plot_persistence_diagram(d)
     plt.show()
 
-- 
cgit v1.2.3


From cb108929433617f553e0a0c7185b3073cce35696 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Thu, 17 Jun 2021 15:43:30 +0200
Subject: Fix #461 and review all error cases (no more prints, warnings and
 exceptions instead)

---
 src/python/CMakeLists.txt                       |   4 +
 src/python/gudhi/persistence_graphical_tools.py | 245 +++++++++++++-----------
 2 files changed, 141 insertions(+), 108 deletions(-)

diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 98f2b85f..1f0d74d4 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -542,6 +542,10 @@ if(PYTHONINTERP_FOUND)
       add_gudhi_py_test(test_dtm_rips_complex)
     endif()
 
+    # persistence graphical tools
+    if(MATPLOTLIB_FOUND)
+      add_gudhi_py_test(test_persistence_graphical_tools)
+    endif()
 
     # Set missing or not modules
     set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES")
diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index b129b375..9bfc19e0 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -13,6 +13,8 @@ from math import isfinite
 import numpy as np
 from functools import lru_cache
 import warnings
+import errno
+import os
 
 from gudhi.reader_utils import read_persistence_intervals_in_dimension
 from gudhi.reader_utils import read_persistence_intervals_grouped_by_dimension
@@ -45,6 +47,9 @@ def __min_birth_max_death(persistence, band=0.0):
             min_birth = float(interval[1][0])
     if band > 0.0:
         max_death += band
+    # can happen if only points at inf death
+    if min_birth == max_death:
+        max_death = max_death + 1.
     return (min_birth, max_death)
 
 
@@ -54,7 +59,7 @@ def _array_handler(a):
                 persistence-compatible list (padding with 0), so that the
                 plot can be performed seamlessly.
     '''
-    if isinstance(a[0][1], np.float64) or isinstance(a[0][1], float):
+    if isinstance(a[0][1], np.floating) or isinstance(a[0][1], float):
         return [[0, x] for x in a]
     else:
         return a
@@ -88,7 +93,7 @@ def _matplotlib_can_use_tex():
         from matplotlib import checkdep_usetex
         return checkdep_usetex(True)
     except ImportError:
-        print("This function is not available, you may be missing matplotlib.")
+        warnings.warn("This function is not available, you may be missing matplotlib.")
 
 
 def plot_persistence_barcode(
@@ -157,53 +162,58 @@ def plot_persistence_barcode(
                     for persistence_interval in diag[key]:
                         persistence.append((key, persistence_interval))
             else:
-                print("file " + persistence_file + " not found.")
-                return None
+                raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file)
 
-        persistence = _array_handler(persistence)
-
-        persistence = _limit_to_max_intervals(persistence, max_intervals,
-                                              key = lambda life_time: life_time[1][1] - life_time[1][0])
-
-        if colormap == None:
-            colormap = plt.cm.Set1.colors
-        if axes == None:
-            _, axes = plt.subplots(1, 1)
-
-        persistence = sorted(persistence, key=lambda birth: birth[1][0])
+        try:
+            persistence = _array_handler(persistence)
+            persistence = _limit_to_max_intervals(persistence, max_intervals,
+                                                  key = lambda life_time: life_time[1][1] - life_time[1][0])
+            (min_birth, max_death) = __min_birth_max_death(persistence)
+            persistence = sorted(persistence, key=lambda birth: birth[1][0])
+        except IndexError:
+            min_birth, max_death = 0., 1.
+            pass
 
-        (min_birth, max_death) = __min_birth_max_death(persistence)
-        ind = 0
         delta = (max_death - min_birth) * inf_delta
         # Replace infinity values with max_death + delta for bar code to be more
         # readable
         infinity = max_death + delta
         axis_start = min_birth - delta
+
+        if axes == None:
+            _, axes = plt.subplots(1, 1)
+        if colormap == None:
+            colormap = plt.cm.Set1.colors
+        ind = 0
+
         # Draw horizontal bars in loop
         for interval in reversed(persistence):
-            if float(interval[1][1]) != float("inf"):
-                # Finite death case
-                axes.barh(
-                    ind,
-                    (interval[1][1] - interval[1][0]),
-                    height=0.8,
-                    left=interval[1][0],
-                    alpha=alpha,
-                    color=colormap[interval[0]],
-                    linewidth=0,
-                )
-            else:
-                # Infinite death case for diagram to be nicer
-                axes.barh(
-                    ind,
-                    (infinity - interval[1][0]),
-                    height=0.8,
-                    left=interval[1][0],
-                    alpha=alpha,
-                    color=colormap[interval[0]],
-                    linewidth=0,
-                )
-            ind = ind + 1
+            try:
+                if float(interval[1][1]) != float("inf"):
+                    # Finite death case
+                    axes.barh(
+                        ind,
+                        (interval[1][1] - interval[1][0]),
+                        height=0.8,
+                        left=interval[1][0],
+                        alpha=alpha,
+                        color=colormap[interval[0]],
+                        linewidth=0,
+                    )
+                else:
+                    # Infinite death case for diagram to be nicer
+                    axes.barh(
+                        ind,
+                        (infinity - interval[1][0]),
+                        height=0.8,
+                        left=interval[1][0],
+                        alpha=alpha,
+                        color=colormap[interval[0]],
+                        linewidth=0,
+                    )
+                ind = ind + 1
+            except IndexError:
+                pass
 
         if legend:
             dimensions = list(set(item[0] for item in persistence))
@@ -218,11 +228,12 @@ def plot_persistence_barcode(
         axes.set_title("Persistence barcode", fontsize=fontsize)
 
         # Ends plot on infinity value and starts a little bit before min_birth
-        axes.axis([axis_start, infinity, 0, ind])
+        if ind != 0:
+            axes.axis([axis_start, infinity, 0, ind])
         return axes
 
     except ImportError:
-        print("This function is not available, you may be missing matplotlib.")
+        warnings.warn("This function is not available, you may be missing matplotlib.")
 
 
 def plot_persistence_diagram(
@@ -296,20 +307,17 @@ def plot_persistence_diagram(
                     for persistence_interval in diag[key]:
                         persistence.append((key, persistence_interval))
             else:
-                print("file " + persistence_file + " not found.")
-                return None
-
-        persistence = _array_handler(persistence)
-
-        persistence = _limit_to_max_intervals(persistence, max_intervals,
-                                              key = lambda life_time: life_time[1][1] - life_time[1][0])
+                raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file)
 
-        if colormap == None:
-            colormap = plt.cm.Set1.colors
-        if axes == None:
-            _, axes = plt.subplots(1, 1)
+        try:
+            persistence = _array_handler(persistence)
+            persistence = _limit_to_max_intervals(persistence, max_intervals,
+                                                  key = lambda life_time: life_time[1][1] - life_time[1][0])
+            min_birth, max_death = __min_birth_max_death(persistence, band)
+        except IndexError:
+            min_birth, max_death = 0., 1.
+            pass
 
-        (min_birth, max_death) = __min_birth_max_death(persistence, band)
         delta = (max_death - min_birth) * inf_delta
         # Replace infinity values with max_death + delta for diagram to be more
         # readable
@@ -317,33 +325,43 @@ def plot_persistence_diagram(
         axis_end = max_death + delta / 2
         axis_start = min_birth - delta
 
+        if axes == None:
+            _, axes = plt.subplots(1, 1)
+        if colormap == None:
+            colormap = plt.cm.Set1.colors
         # bootstrap band
         if band > 0.0:
             x = np.linspace(axis_start, infinity, 1000)
             axes.fill_between(x, x, x + band, alpha=alpha, facecolor="red")
         # lower diag patch
         if greyblock:
-            axes.add_patch(mpatches.Polygon([[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]], fill=True, color='lightgrey'))
+            axes.add_patch(mpatches.Polygon([[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]],
+                                            fill=True, color='lightgrey'))
+        # line display of equation : birth = death
+        axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
+
         # Draw points in loop
         pts_at_infty = False  # Records presence of pts at infty
         for interval in reversed(persistence):
-            if float(interval[1][1]) != float("inf"):
-                # Finite death case
-                axes.scatter(
-                    interval[1][0],
-                    interval[1][1],
-                    alpha=alpha,
-                    color=colormap[interval[0]],
-                )
-            else:
-                pts_at_infty = True
-                # Infinite death case for diagram to be nicer
-                axes.scatter(
-                    interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]]
-                )
+            try:
+                if float(interval[1][1]) != float("inf"):
+                    # Finite death case
+                    axes.scatter(
+                        interval[1][0],
+                        interval[1][1],
+                        alpha=alpha,
+                        color=colormap[interval[0]],
+                    )
+                else:
+                    pts_at_infty = True
+                    # Infinite death case for diagram to be nicer
+                    axes.scatter(
+                        interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]]
+                    )
+            except IndexError:
+                pass
         if pts_at_infty:
             # infinity line and text
-            axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
             axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha)
             # Infinity label
             yt = axes.get_yticks()
@@ -371,7 +389,7 @@ def plot_persistence_diagram(
         return axes
 
     except ImportError:
-        print("This function is not available, you may be missing matplotlib.")
+        warnings.warn("This function is not available, you may be missing matplotlib.")
 
 
 def plot_persistence_density(
@@ -461,27 +479,7 @@ def plot_persistence_density(
                     persistence_file=persistence_file, only_this_dim=dimension
                 )
             else:
-                print("file " + persistence_file + " not found.")
-                return None
-
-        if len(persistence) > 0:
-            persistence = _array_handler(persistence)
-            persistence_dim = np.array(
-                [
-                    (dim_interval[1][0], dim_interval[1][1])
-                    for dim_interval in persistence
-                    if (dim_interval[0] == dimension) or (dimension is None)
-                ]
-            )
-
-        persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])]
-
-        persistence_dim = np.array(_limit_to_max_intervals(persistence_dim, max_intervals,
-                                                           key = lambda life_time: life_time[1] - life_time[0]))
-
-        # Set as numpy array birth and death (remove undefined values - inf and NaN)
-        birth = persistence_dim[:, 0]
-        death = persistence_dim[:, 1]
+                raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file)
 
         # default cmap value cannot be done at argument definition level as matplotlib is not yet defined.
         if cmap is None:
@@ -489,23 +487,56 @@ def plot_persistence_density(
         if axes == None:
             _, axes = plt.subplots(1, 1)
 
+        try:
+            if len(persistence) > 0:
+                # if not read from file but given by an argument
+                persistence = _array_handler(persistence)
+                persistence_dim = np.array(
+                    [
+                        (dim_interval[1][0], dim_interval[1][1])
+                        for dim_interval in persistence
+                        if (dim_interval[0] == dimension) or (dimension is None)
+                    ]
+                )
+            persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])]
+            persistence_dim = np.array(_limit_to_max_intervals(persistence_dim, max_intervals,
+                                                               key = lambda life_time: life_time[1] - life_time[0]))
+
+            # Set as numpy array birth and death (remove undefined values - inf and NaN)
+            birth = persistence_dim[:, 0]
+            death = persistence_dim[:, 1]
+            birth_min = birth.min()
+            birth_max = birth.max()
+            death_min = death.min()
+            death_max = death.max()
+
+            # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
+            k = kde.gaussian_kde([birth, death], bw_method=bw_method)
+            xi, yi = np.mgrid[
+                birth_min : birth_max : nbins * 1j,
+                death_min : death_max : nbins * 1j,
+            ]
+            zi = k(np.vstack([xi.flatten(), yi.flatten()]))
+            # Make the plot
+            img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading='auto')
+
+        # IndexError on empty diagrams, ValueError on only inf death values
+        except (IndexError, ValueError):
+            birth_min = 0.
+            birth_max = 1.
+            death_min = 0.
+            death_max = 1.
+            pass
+
         # line display of equation : birth = death
-        x = np.linspace(death.min(), birth.max(), 1000)
+        x = np.linspace(death_min, birth_max, 1000)
         axes.plot(x, x, color="k", linewidth=1.0)
 
-        # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
-        k = kde.gaussian_kde([birth, death], bw_method=bw_method)
-        xi, yi = np.mgrid[
-            birth.min() : birth.max() : nbins * 1j,
-            death.min() : death.max() : nbins * 1j,
-        ]
-        zi = k(np.vstack([xi.flatten(), yi.flatten()]))
-
-        # Make the plot
-        img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading='auto')
-
         if greyblock:
-            axes.add_patch(mpatches.Polygon([[birth.min(), birth.min()], [death.max(), birth.min()], [death.max(), death.max()]], fill=True, color='lightgrey'))
+            axes.add_patch(mpatches.Polygon([[birth_min, birth_min],
+                                             [death_max, birth_min],
+                                             [death_max, death_max]],
+                                            fill=True, color='lightgrey'))
 
         if legend:
             plt.colorbar(img, ax=axes)
@@ -517,6 +548,4 @@ def plot_persistence_density(
         return axes
 
     except ImportError:
-        print(
-            "This function is not available, you may be missing matplotlib and/or scipy."
-        )
+        warnings.warn("This function is not available, you may be missing matplotlib and/or scipy.")
-- 
cgit v1.2.3


From 72f72770ccf0d1ddb78a7f23103b2777f407e72c Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Thu, 17 Jun 2021 16:29:18 +0200
Subject: Forgot to add test file for graphical

---
 .../test/test_persistence_graphical_tools.py       | 105 +++++++++++++++++++++
 1 file changed, 105 insertions(+)
 create mode 100644 src/python/test/test_persistence_graphical_tools.py

diff --git a/src/python/test/test_persistence_graphical_tools.py b/src/python/test/test_persistence_graphical_tools.py
new file mode 100644
index 00000000..1ad1ae23
--- /dev/null
+++ b/src/python/test/test_persistence_graphical_tools.py
@@ -0,0 +1,105 @@
+""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+    Author(s):       Vincent Rouvreau
+
+    Copyright (C) 2021 Inria
+
+    Modification(s):
+      - YYYY/MM Author: Description of the modification
+"""
+
+import gudhi as gd
+import numpy as np
+import matplotlib as plt
+import pytest
+
+def test_array_handler():
+    diags = np.array([[1, 2], [3, 4], [5, 6]], np.float)
+    arr_diags = gd.persistence_graphical_tools._array_handler(diags)
+    for idx in range(len(diags)):
+        assert arr_diags[idx][0] == 0
+        np.testing.assert_array_equal(arr_diags[idx][1], diags[idx])
+
+    diags = [(1., 2.), (3., 4.), (5., 6.)]
+    arr_diags = gd.persistence_graphical_tools._array_handler(diags)
+    for idx in range(len(diags)):
+        assert arr_diags[idx][0] == 0
+        assert arr_diags[idx][1] == diags[idx]
+
+    diags = [(0, (1., 2.)), (0, (3., 4.)), (0, (5., 6.))]
+    assert gd.persistence_graphical_tools._array_handler(diags) == diags
+
+def test_min_birth_max_death():
+    diags = [
+        (0, (0., float("inf"))),
+        (0, (0.0983494, float("inf"))),
+        (0, (0., 0.122545)),
+        (0, (0., 0.12047)),
+        (0, (0., 0.118398)),
+        (0, (0.118398, 1.)),
+        (0, (0., 0.117908)),
+        (0, (0., 0.112307)),
+        (0, (0., 0.107535)),
+        (0, (0., 0.106382)),
+    ]
+    assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (0., 1.)
+    assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band=4.) == (0., 5.)
+
+def test_limit_min_birth_max_death():
+    diags = [
+        (0, (2., float("inf"))),
+        (0, (2., float("inf"))),
+    ]
+    assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (2., 3.)
+    assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band = 4.) == (2., 6.)
+
+def test_limit_to_max_intervals():
+    diags = [
+        (0, (0., float("inf"))),
+        (0, (0.0983494, float("inf"))),
+        (0, (0., 0.122545)),
+        (0, (0., 0.12047)),
+        (0, (0., 0.118398)),
+        (0, (0.118398, 1.)),
+        (0, (0., 0.117908)),
+        (0, (0., 0.112307)),
+        (0, (0., 0.107535)),
+        (0, (0., 0.106382)),
+    ]
+    # check no warnings if max_intervals equals to the diagrams number
+    with pytest.warns(None) as record:
+        truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals(diags, 10,
+            key = lambda life_time: life_time[1][1] - life_time[1][0])
+        # check diagrams are not sorted
+        assert truncated_diags == diags
+    assert len(record) == 0
+
+    # check warning if max_intervals lower than the diagrams number
+    with pytest.warns(UserWarning) as record:
+        truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals(diags, 5,
+            key = lambda life_time: life_time[1][1] - life_time[1][0])
+        # check diagrams are truncated and sorted by life time
+        assert truncated_diags == [(0, (0.,         float("inf"))),
+                                   (0, (0.0983494, float("inf"))),
+                                   (0, (0.118398,  1.0)),
+                                   (0, (0.,         0.122545)),
+                                   (0, (0.,         0.12047))]
+    assert len(record) == 1
+
+def _limit_plot_persistence(function):
+    pplot = function(persistence=[()])
+    assert issubclass(type(pplot), plt.axes.SubplotBase)
+    pplot = function(persistence=[(0, float("inf"))])
+    assert issubclass(type(pplot), plt.axes.SubplotBase)
+
+def test_limit_plot_persistence():
+    for function in [gd.plot_persistence_barcode, gd.plot_persistence_diagram, gd.plot_persistence_density]:
+        _limit_plot_persistence(function)
+
+def _non_existing_persistence_file(function):
+    with pytest.raises(FileNotFoundError):
+        function(persistence_file="pouetpouettralala.toubiloubabdou")
+
+def test_non_existing_persistence_file():
+    for function in [gd.plot_persistence_barcode, gd.plot_persistence_diagram, gd.plot_persistence_density]:
+        _non_existing_persistence_file(function)
-- 
cgit v1.2.3


From f9b1e50a6adaadf88c4940cb4214f9ecef542144 Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 22 Jun 2021 18:47:46 +0200
Subject: black -l 120 modified files

---
 src/python/gudhi/persistence_graphical_tools.py    | 144 +++++++++++----------
 .../test/test_persistence_graphical_tools.py       |  82 +++++++-----
 2 files changed, 120 insertions(+), 106 deletions(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 9bfc19e0..837218ca 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -25,6 +25,7 @@ __license__ = "MIT"
 
 _gudhi_matplotlib_use_tex = True
 
+
 def __min_birth_max_death(persistence, band=0.0):
     """This function returns (min_birth, max_death) from the persistence.
 
@@ -49,23 +50,24 @@ def __min_birth_max_death(persistence, band=0.0):
         max_death += band
     # can happen if only points at inf death
     if min_birth == max_death:
-        max_death = max_death + 1.
+        max_death = max_death + 1.0
     return (min_birth, max_death)
 
 
 def _array_handler(a):
-    '''
+    """
     :param a: if array, assumes it is a (n x 2) np.array and return a
                 persistence-compatible list (padding with 0), so that the
                 plot can be performed seamlessly.
-    '''
+    """
     if isinstance(a[0][1], np.floating) or isinstance(a[0][1], float):
         return [[0, x] for x in a]
     else:
         return a
 
+
 def _limit_to_max_intervals(persistence, max_intervals, key):
-    '''This function returns truncated persistence if length is bigger than max_intervals.
+    """This function returns truncated persistence if length is bigger than max_intervals.
     :param persistence: Persistence intervals values list. Can be grouped by dimension or not.
     :type persistence: an array of (dimension, array of (birth, death)) or an array of (birth, death).
     :param max_intervals: maximal number of intervals to display.
@@ -74,16 +76,18 @@ def _limit_to_max_intervals(persistence, max_intervals, key):
     :type max_intervals: int.
     :param key: key function for sort algorithm.
     :type key: function or lambda.
-    '''
+    """
     if max_intervals > 0 and max_intervals < len(persistence):
-        warnings.warn('There are %s intervals given as input, whereas max_intervals is set to %s.' %
-                        (len(persistence), max_intervals)
-                     )
+        warnings.warn(
+            "There are %s intervals given as input, whereas max_intervals is set to %s."
+            % (len(persistence), max_intervals)
+        )
         # Sort by life time, then takes only the max_intervals elements
-        return sorted(persistence, key = key, reverse = True)[:max_intervals]
+        return sorted(persistence, key=key, reverse=True)[:max_intervals]
     else:
         return persistence
 
+
 @lru_cache(maxsize=1)
 def _matplotlib_can_use_tex():
     """This function returns True if matplotlib can deal with LaTeX, False otherwise.
@@ -91,6 +95,7 @@ def _matplotlib_can_use_tex():
     """
     try:
         from matplotlib import checkdep_usetex
+
         return checkdep_usetex(True)
     except ImportError:
         warnings.warn("This function is not available, you may be missing matplotlib.")
@@ -144,20 +149,19 @@ def plot_persistence_barcode(
         import matplotlib.pyplot as plt
         import matplotlib.patches as mpatches
         from matplotlib import rc
+
         if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex():
-            plt.rc('text', usetex=True)
-            plt.rc('font', family='serif')
+            plt.rc("text", usetex=True)
+            plt.rc("font", family="serif")
         else:
-            plt.rc('text', usetex=False)
-            plt.rc('font', family='DejaVu Sans')
+            plt.rc("text", usetex=False)
+            plt.rc("font", family="DejaVu Sans")
 
         if persistence_file != "":
             if path.isfile(persistence_file):
                 # Reset persistence
                 persistence = []
-                diag = read_persistence_intervals_grouped_by_dimension(
-                    persistence_file=persistence_file
-                )
+                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))
@@ -166,12 +170,13 @@ def plot_persistence_barcode(
 
         try:
             persistence = _array_handler(persistence)
-            persistence = _limit_to_max_intervals(persistence, max_intervals,
-                                                  key = lambda life_time: life_time[1][1] - life_time[1][0])
+            persistence = _limit_to_max_intervals(
+                persistence, max_intervals, key=lambda life_time: life_time[1][1] - life_time[1][0]
+            )
             (min_birth, max_death) = __min_birth_max_death(persistence)
             persistence = sorted(persistence, key=lambda birth: birth[1][0])
         except IndexError:
-            min_birth, max_death = 0., 1.
+            min_birth, max_death = 0.0, 1.0
             pass
 
         delta = (max_death - min_birth) * inf_delta
@@ -218,11 +223,7 @@ def plot_persistence_barcode(
         if legend:
             dimensions = list(set(item[0] for item in persistence))
             axes.legend(
-                handles=[
-                    mpatches.Patch(color=colormap[dim], label=str(dim))
-                    for dim in dimensions
-                ],
-                loc="lower right",
+                handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions], loc="lower right",
             )
 
         axes.set_title("Persistence barcode", fontsize=fontsize)
@@ -247,7 +248,7 @@ def plot_persistence_diagram(
     colormap=None,
     axes=None,
     fontsize=16,
-    greyblock=True
+    greyblock=True,
 ):
     """This function plots the persistence diagram from persistence values
     list, a np.array of shape (N x 2) representing a diagram in a single
@@ -289,20 +290,19 @@ def plot_persistence_diagram(
         import matplotlib.pyplot as plt
         import matplotlib.patches as mpatches
         from matplotlib import rc
+
         if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex():
-            plt.rc('text', usetex=True)
-            plt.rc('font', family='serif')
+            plt.rc("text", usetex=True)
+            plt.rc("font", family="serif")
         else:
-            plt.rc('text', usetex=False)
-            plt.rc('font', family='DejaVu Sans')
+            plt.rc("text", usetex=False)
+            plt.rc("font", family="DejaVu Sans")
 
         if persistence_file != "":
             if path.isfile(persistence_file):
                 # Reset persistence
                 persistence = []
-                diag = read_persistence_intervals_grouped_by_dimension(
-                    persistence_file=persistence_file
-                )
+                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))
@@ -311,11 +311,12 @@ def plot_persistence_diagram(
 
         try:
             persistence = _array_handler(persistence)
-            persistence = _limit_to_max_intervals(persistence, max_intervals,
-                                                  key = lambda life_time: life_time[1][1] - life_time[1][0])
+            persistence = _limit_to_max_intervals(
+                persistence, max_intervals, key=lambda life_time: life_time[1][1] - life_time[1][0]
+            )
             min_birth, max_death = __min_birth_max_death(persistence, band)
         except IndexError:
-            min_birth, max_death = 0., 1.
+            min_birth, max_death = 0.0, 1.0
             pass
 
         delta = (max_death - min_birth) * inf_delta
@@ -335,8 +336,13 @@ def plot_persistence_diagram(
             axes.fill_between(x, x, x + band, alpha=alpha, facecolor="red")
         # lower diag patch
         if greyblock:
-            axes.add_patch(mpatches.Polygon([[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]],
-                                            fill=True, color='lightgrey'))
+            axes.add_patch(
+                mpatches.Polygon(
+                    [[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]],
+                    fill=True,
+                    color="lightgrey",
+                )
+            )
         # line display of equation : birth = death
         axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
 
@@ -347,17 +353,12 @@ def plot_persistence_diagram(
                 if float(interval[1][1]) != float("inf"):
                     # Finite death case
                     axes.scatter(
-                        interval[1][0],
-                        interval[1][1],
-                        alpha=alpha,
-                        color=colormap[interval[0]],
+                        interval[1][0], interval[1][1], alpha=alpha, color=colormap[interval[0]],
                     )
                 else:
                     pts_at_infty = True
                     # Infinite death case for diagram to be nicer
-                    axes.scatter(
-                        interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]]
-                    )
+                    axes.scatter(interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]])
             except IndexError:
                 pass
         if pts_at_infty:
@@ -365,27 +366,22 @@ def plot_persistence_diagram(
             axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha)
             # Infinity label
             yt = axes.get_yticks()
-            yt = yt[np.where(yt < axis_end)] # to avoid ploting ticklabel higher than infinity
+            yt = yt[np.where(yt < axis_end)]  # to avoid ploting ticklabel higher than infinity
             yt = np.append(yt, infinity)
             ytl = ["%.3f" % e for e in yt]  # to avoid float precision error
-            ytl[-1] = r'$+\infty$'
+            ytl[-1] = r"$+\infty$"
             axes.set_yticks(yt)
             axes.set_yticklabels(ytl)
 
         if legend:
             dimensions = list(set(item[0] for item in persistence))
-            axes.legend(
-                handles=[
-                    mpatches.Patch(color=colormap[dim], label=str(dim))
-                    for dim in dimensions
-                ]
-            )
+            axes.legend(handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions])
 
         axes.set_xlabel("Birth", fontsize=fontsize)
         axes.set_ylabel("Death", fontsize=fontsize)
         axes.set_title("Persistence diagram", fontsize=fontsize)
         # Ends plot on infinity value and starts a little bit before min_birth
-        axes.axis([axis_start, axis_end, axis_start, infinity + delta/2])
+        axes.axis([axis_start, axis_end, axis_start, infinity + delta / 2])
         return axes
 
     except ImportError:
@@ -403,7 +399,7 @@ def plot_persistence_density(
     legend=False,
     axes=None,
     fontsize=16,
-    greyblock=False
+    greyblock=False,
 ):
     """This function plots the persistence density from persistence
     values list, np.array of shape (N x 2) representing a diagram
@@ -463,12 +459,13 @@ def plot_persistence_density(
         import matplotlib.patches as mpatches
         from scipy.stats import kde
         from matplotlib import rc
+
         if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex():
-            plt.rc('text', usetex=True)
-            plt.rc('font', family='serif')
+            plt.rc("text", usetex=True)
+            plt.rc("font", family="serif")
         else:
-            plt.rc('text', usetex=False)
-            plt.rc('font', family='DejaVu Sans')
+            plt.rc("text", usetex=False)
+            plt.rc("font", family="DejaVu Sans")
 
         if persistence_file != "":
             if dimension is None:
@@ -499,8 +496,11 @@ def plot_persistence_density(
                     ]
                 )
             persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])]
-            persistence_dim = np.array(_limit_to_max_intervals(persistence_dim, max_intervals,
-                                                               key = lambda life_time: life_time[1] - life_time[0]))
+            persistence_dim = np.array(
+                _limit_to_max_intervals(
+                    persistence_dim, max_intervals, key=lambda life_time: life_time[1] - life_time[0]
+                )
+            )
 
             # Set as numpy array birth and death (remove undefined values - inf and NaN)
             birth = persistence_dim[:, 0]
@@ -513,19 +513,18 @@ def plot_persistence_density(
             # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
             k = kde.gaussian_kde([birth, death], bw_method=bw_method)
             xi, yi = np.mgrid[
-                birth_min : birth_max : nbins * 1j,
-                death_min : death_max : nbins * 1j,
+                birth_min : birth_max : nbins * 1j, death_min : death_max : nbins * 1j,
             ]
             zi = k(np.vstack([xi.flatten(), yi.flatten()]))
             # Make the plot
-            img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading='auto')
+            img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading="auto")
 
         # IndexError on empty diagrams, ValueError on only inf death values
         except (IndexError, ValueError):
-            birth_min = 0.
-            birth_max = 1.
-            death_min = 0.
-            death_max = 1.
+            birth_min = 0.0
+            birth_max = 1.0
+            death_min = 0.0
+            death_max = 1.0
             pass
 
         # line display of equation : birth = death
@@ -533,10 +532,13 @@ def plot_persistence_density(
         axes.plot(x, x, color="k", linewidth=1.0)
 
         if greyblock:
-            axes.add_patch(mpatches.Polygon([[birth_min, birth_min],
-                                             [death_max, birth_min],
-                                             [death_max, death_max]],
-                                            fill=True, color='lightgrey'))
+            axes.add_patch(
+                mpatches.Polygon(
+                    [[birth_min, birth_min], [death_max, birth_min], [death_max, death_max]],
+                    fill=True,
+                    color="lightgrey",
+                )
+            )
 
         if legend:
             plt.colorbar(img, ax=axes)
diff --git a/src/python/test/test_persistence_graphical_tools.py b/src/python/test/test_persistence_graphical_tools.py
index 1ad1ae23..7d9bae90 100644
--- a/src/python/test/test_persistence_graphical_tools.py
+++ b/src/python/test/test_persistence_graphical_tools.py
@@ -13,6 +13,7 @@ import numpy as np
 import matplotlib as plt
 import pytest
 
+
 def test_array_handler():
     diags = np.array([[1, 2], [3, 4], [5, 6]], np.float)
     arr_diags = gd.persistence_graphical_tools._array_handler(diags)
@@ -20,86 +21,97 @@ def test_array_handler():
         assert arr_diags[idx][0] == 0
         np.testing.assert_array_equal(arr_diags[idx][1], diags[idx])
 
-    diags = [(1., 2.), (3., 4.), (5., 6.)]
+    diags = [(1.0, 2.0), (3.0, 4.0), (5.0, 6.0)]
     arr_diags = gd.persistence_graphical_tools._array_handler(diags)
     for idx in range(len(diags)):
         assert arr_diags[idx][0] == 0
         assert arr_diags[idx][1] == diags[idx]
 
-    diags = [(0, (1., 2.)), (0, (3., 4.)), (0, (5., 6.))]
+    diags = [(0, (1.0, 2.0)), (0, (3.0, 4.0)), (0, (5.0, 6.0))]
     assert gd.persistence_graphical_tools._array_handler(diags) == diags
 
+
 def test_min_birth_max_death():
     diags = [
-        (0, (0., float("inf"))),
+        (0, (0.0, float("inf"))),
         (0, (0.0983494, float("inf"))),
-        (0, (0., 0.122545)),
-        (0, (0., 0.12047)),
-        (0, (0., 0.118398)),
-        (0, (0.118398, 1.)),
-        (0, (0., 0.117908)),
-        (0, (0., 0.112307)),
-        (0, (0., 0.107535)),
-        (0, (0., 0.106382)),
+        (0, (0.0, 0.122545)),
+        (0, (0.0, 0.12047)),
+        (0, (0.0, 0.118398)),
+        (0, (0.118398, 1.0)),
+        (0, (0.0, 0.117908)),
+        (0, (0.0, 0.112307)),
+        (0, (0.0, 0.107535)),
+        (0, (0.0, 0.106382)),
     ]
-    assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (0., 1.)
-    assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band=4.) == (0., 5.)
+    assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (0.0, 1.0)
+    assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band=4.0) == (0.0, 5.0)
+
 
 def test_limit_min_birth_max_death():
     diags = [
-        (0, (2., float("inf"))),
-        (0, (2., float("inf"))),
+        (0, (2.0, float("inf"))),
+        (0, (2.0, float("inf"))),
     ]
-    assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (2., 3.)
-    assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band = 4.) == (2., 6.)
+    assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (2.0, 3.0)
+    assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band=4.0) == (2.0, 6.0)
+
 
 def test_limit_to_max_intervals():
     diags = [
-        (0, (0., float("inf"))),
+        (0, (0.0, float("inf"))),
         (0, (0.0983494, float("inf"))),
-        (0, (0., 0.122545)),
-        (0, (0., 0.12047)),
-        (0, (0., 0.118398)),
-        (0, (0.118398, 1.)),
-        (0, (0., 0.117908)),
-        (0, (0., 0.112307)),
-        (0, (0., 0.107535)),
-        (0, (0., 0.106382)),
+        (0, (0.0, 0.122545)),
+        (0, (0.0, 0.12047)),
+        (0, (0.0, 0.118398)),
+        (0, (0.118398, 1.0)),
+        (0, (0.0, 0.117908)),
+        (0, (0.0, 0.112307)),
+        (0, (0.0, 0.107535)),
+        (0, (0.0, 0.106382)),
     ]
     # check no warnings if max_intervals equals to the diagrams number
     with pytest.warns(None) as record:
-        truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals(diags, 10,
-            key = lambda life_time: life_time[1][1] - life_time[1][0])
+        truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals(
+            diags, 10, key=lambda life_time: life_time[1][1] - life_time[1][0]
+        )
         # check diagrams are not sorted
         assert truncated_diags == diags
     assert len(record) == 0
 
     # check warning if max_intervals lower than the diagrams number
     with pytest.warns(UserWarning) as record:
-        truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals(diags, 5,
-            key = lambda life_time: life_time[1][1] - life_time[1][0])
+        truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals(
+            diags, 5, key=lambda life_time: life_time[1][1] - life_time[1][0]
+        )
         # check diagrams are truncated and sorted by life time
-        assert truncated_diags == [(0, (0.,         float("inf"))),
-                                   (0, (0.0983494, float("inf"))),
-                                   (0, (0.118398,  1.0)),
-                                   (0, (0.,         0.122545)),
-                                   (0, (0.,         0.12047))]
+        assert truncated_diags == [
+            (0, (0.0, float("inf"))),
+            (0, (0.0983494, float("inf"))),
+            (0, (0.118398, 1.0)),
+            (0, (0.0, 0.122545)),
+            (0, (0.0, 0.12047)),
+        ]
     assert len(record) == 1
 
+
 def _limit_plot_persistence(function):
     pplot = function(persistence=[()])
     assert issubclass(type(pplot), plt.axes.SubplotBase)
     pplot = function(persistence=[(0, float("inf"))])
     assert issubclass(type(pplot), plt.axes.SubplotBase)
 
+
 def test_limit_plot_persistence():
     for function in [gd.plot_persistence_barcode, gd.plot_persistence_diagram, gd.plot_persistence_density]:
         _limit_plot_persistence(function)
 
+
 def _non_existing_persistence_file(function):
     with pytest.raises(FileNotFoundError):
         function(persistence_file="pouetpouettralala.toubiloubabdou")
 
+
 def test_non_existing_persistence_file():
     for function in [gd.plot_persistence_barcode, gd.plot_persistence_diagram, gd.plot_persistence_density]:
         _non_existing_persistence_file(function)
-- 
cgit v1.2.3


From 87da488d24c70cbd470ad1c2dae762af68cd227e Mon Sep 17 00:00:00 2001
From: ROUVREAU Vincent <vincent.rouvreau@inria.fr>
Date: Tue, 24 Aug 2021 14:50:12 +0200
Subject: code review: isinstance can take several instances types as a list

---
 src/python/gudhi/persistence_graphical_tools.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 837218ca..70d550ab 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -60,7 +60,7 @@ def _array_handler(a):
                 persistence-compatible list (padding with 0), so that the
                 plot can be performed seamlessly.
     """
-    if isinstance(a[0][1], np.floating) or isinstance(a[0][1], float):
+    if isinstance(a[0][1], (np.floating, float)):
         return [[0, x] for x in a]
     else:
         return a
-- 
cgit v1.2.3


From 37a141533397568e7070c734e21ef9c4dc85d132 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Thu, 10 Feb 2022 16:16:14 +0100
Subject: Add SimplexTree copy method and its test

---
 src/python/gudhi/simplex_tree.pxd    |  1 +
 src/python/gudhi/simplex_tree.pyx    | 14 ++++++++++++++
 src/python/test/test_simplex_tree.py | 19 ++++++++++++++++++-
 3 files changed, 33 insertions(+), 1 deletion(-)

diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 006a24ed..92139db4 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -45,6 +45,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
 
     cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>":
         Simplex_tree_interface_full_featured() nogil
+        Simplex_tree_interface_full_featured(Simplex_tree_interface_full_featured&) nogil
         double simplex_filtration(vector[int] simplex) nogil
         void assign_simplex_filtration(vector[int] simplex, double filtration) nogil
         void initialize_filtration() nogil
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index c3720936..6b3116a4 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -63,6 +63,20 @@ cdef class SimplexTree:
          """
         return self.pcohptr != NULL
 
+    def copy(self):
+        """ 
+        :returns: A simplex tree that is a deep copy itself.
+        :rtype: SimplexTree
+        """
+        stree = SimplexTree()
+        cdef Simplex_tree_interface_full_featured* stree_ptr
+        cdef Simplex_tree_interface_full_featured* self_ptr=self.get_ptr()
+        with nogil:
+            stree_ptr = new Simplex_tree_interface_full_featured(dereference(self_ptr))
+
+        stree.thisptr = <intptr_t>(stree_ptr)
+        return stree
+
     def filtration(self, simplex):
         """This function returns the filtration value for a given N-simplex in
         this simplicial complex, or +infinity if it is not in the complex.
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 31c46213..dac45288 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -447,4 +447,21 @@ def test_persistence_intervals_in_dimension():
     assert np.array_equal(H2, np.array([[ 0., float("inf")]]))
     # Test empty case
     assert st.persistence_intervals_in_dimension(3).shape == (0, 2)
-    
\ No newline at end of file
+
+def test_simplex_tree_copy():
+    st = SimplexTree()
+    st .insert([1,2,3], 0.)
+    a = st.copy()
+    # TODO(VR): when #463 is merged, replace with
+    # assert a == st
+    assert a.num_vertices() == st.num_vertices()
+    assert a.num_simplices() == st.num_simplices()
+    st_filt_list = list(st.get_filtration())
+    assert list(a.get_filtration()) == st_filt_list
+
+    a.remove_maximal_simplex([1, 2, 3])
+    a_filt_list = list(a.get_filtration())
+    assert len(a_filt_list) < len(st_filt_list)
+
+    for a_splx in a_filt_list:
+        assert a_splx in st_filt_list
-- 
cgit v1.2.3


From fb8ce008feadcaf6a936740a3ed54d50970c731c Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Fri, 11 Feb 2022 23:11:26 +0100
Subject: __copy__, __deepcopy__, copy, and copy ctors. Still pb with the doc

---
 .github/next_release.md              |   3 ++
 src/python/gudhi/simplex_tree.pyx    |  54 +++++++++++++++----
 src/python/test/test_simplex_tree.py | 102 +++++++++++++++++++++++++++++++----
 3 files changed, 140 insertions(+), 19 deletions(-)

diff --git a/.github/next_release.md b/.github/next_release.md
index e21b25c7..3946404b 100644
--- a/.github/next_release.md
+++ b/.github/next_release.md
@@ -13,6 +13,9 @@ Below is a list of changes made since GUDHI 3.5.0:
 - [Representations](https://gudhi.inria.fr/python/latest/representations.html#gudhi.representations.vector_methods.BettiCurve)
      - A more flexible Betti curve class capable of computing exact curves
 
+- [Simplex tree](https://gudhi.inria.fr/python/latest/simplex_tree_ref.html)
+     - `__copy__`, `__deepcopy__`, `copy` and copy constructors
+
 - Installation
      - Boost &ge; 1.66.0 is now required (was &ge; 1.56.0).
 
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 6b3116a4..ed7c3b92 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -30,6 +30,7 @@ cdef class SimplexTree:
     # unfortunately 'cdef public Simplex_tree_interface_full_featured* thisptr' is not possible
     # Use intptr_t instead to cast the pointer
     cdef public intptr_t thisptr
+    cdef bool __thisptr_to_be_deleted
 
     # Get the pointer casted as it should be
     cdef Simplex_tree_interface_full_featured* get_ptr(self) nogil:
@@ -38,17 +39,36 @@ cdef class SimplexTree:
     cdef Simplex_tree_persistence_interface * pcohptr
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self):
+    def __init__(self, other = None, copy = True):
         """SimplexTree constructor.
+        :param other: If `other` is a SimplexTree (default = None), the SimplexTree is constructed from a deep/shallow copy of `other`.
+        :type other: SimplexTree
+        :param copy: If `True`, the copy will be deep and if `False, the copy will be shallow. Default is `True`.
+        :type copy: bool
+        :returns: A simplex tree that is a (deep or shallow) copy of itself.
+        :rtype: SimplexTree
+        :note: copy constructor requires :func:`compute_persistence` to be launched again as the result is not copied.
         """
 
     # The real cython constructor
-    def __cinit__(self):
-        self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured())
+    def __cinit__(self, other = None, copy = True):
+        cdef SimplexTree ostr
+        if other and type(other) is SimplexTree:
+            ostr = <SimplexTree> other
+            if copy:
+                self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured(dereference(ostr.get_ptr())))
+            else:
+                self.thisptr = ostr.thisptr
+                # Avoid double free - The original is in charge of deletion
+                self.__thisptr_to_be_deleted = False
+        else:
+            self.__thisptr_to_be_deleted = True
+            self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured())
 
     def __dealloc__(self):
         cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr()
-        if ptr != NULL:
+        # Avoid double free - The original is in charge of deletion
+        if ptr != NULL and self.__thisptr_to_be_deleted:
             del ptr
         if self.pcohptr != NULL:
             del self.pcohptr
@@ -63,20 +83,34 @@ cdef class SimplexTree:
          """
         return self.pcohptr != NULL
 
-    def copy(self):
+    def copy(self, deep=True):
         """ 
-        :returns: A simplex tree that is a deep copy itself.
+        :param deep: If `True`, the copy will be deep and if `False`, the copy will be shallow. Default is `True`.
+        :type deep: bool
+        :returns: A simplex tree that is a (deep or shallow) copy of itself.
         :rtype: SimplexTree
+        :note: copy requires :func:`compute_persistence` to be launched again as the result is not copied.
         """
         stree = SimplexTree()
         cdef Simplex_tree_interface_full_featured* stree_ptr
         cdef Simplex_tree_interface_full_featured* self_ptr=self.get_ptr()
-        with nogil:
-            stree_ptr = new Simplex_tree_interface_full_featured(dereference(self_ptr))
-
-        stree.thisptr = <intptr_t>(stree_ptr)
+        if deep:
+            with nogil:
+                stree_ptr = new Simplex_tree_interface_full_featured(dereference(self_ptr))
+    
+            stree.thisptr = <intptr_t>(stree_ptr)
+        else:
+            stree.thisptr = self.thisptr
+            # Avoid double free - The original is in charge of deletion
+            stree.__thisptr_to_be_deleted = False
         return stree
 
+    def __copy__(self):
+        return self.copy(deep=False)
+
+    def __deepcopy__(self):
+        return self.copy(deep=True)
+
     def filtration(self, simplex):
         """This function returns the filtration value for a given N-simplex in
         this simplicial complex, or +infinity if it is not in the complex.
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index dac45288..6db6d8fb 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -448,20 +448,104 @@ def test_persistence_intervals_in_dimension():
     # Test empty case
     assert st.persistence_intervals_in_dimension(3).shape == (0, 2)
 
-def test_simplex_tree_copy():
+def test_simplex_tree_deep_copy():
     st = SimplexTree()
-    st .insert([1,2,3], 0.)
-    a = st.copy()
+    st.insert([1, 2, 3], 0.)
+    # persistence is not copied
+    st.compute_persistence()
+
+    st_copy = st.copy(deep=True)
     # TODO(VR): when #463 is merged, replace with
-    # assert a == st
-    assert a.num_vertices() == st.num_vertices()
-    assert a.num_simplices() == st.num_simplices()
+    # assert st_copy == st
+    assert st_copy.num_vertices() == st.num_vertices()
+    assert st_copy.num_simplices() == st.num_simplices()
     st_filt_list = list(st.get_filtration())
-    assert list(a.get_filtration()) == st_filt_list
+    assert list(st_copy.get_filtration()) == st_filt_list
+
+    assert st.__is_persistence_defined() == True
+    assert st_copy.__is_persistence_defined() == False
 
-    a.remove_maximal_simplex([1, 2, 3])
-    a_filt_list = list(a.get_filtration())
+    st_copy.remove_maximal_simplex([1, 2, 3])
+    a_filt_list = list(st_copy.get_filtration())
     assert len(a_filt_list) < len(st_filt_list)
 
     for a_splx in a_filt_list:
         assert a_splx in st_filt_list
+    
+    # test double free
+    del st
+    del st_copy
+
+def test_simplex_tree_shallow_copy():
+    st = SimplexTree()
+    st.insert([1, 2, 3], 0.)
+    # persistence is not copied
+    st.compute_persistence()
+
+    st_copy = st.copy(deep=False)
+    # TODO(VR): when #463 is merged, replace with
+    # assert st_copy == st
+    assert st_copy.num_vertices() == st.num_vertices()
+    assert st_copy.num_simplices() == st.num_simplices()
+    assert list(st_copy.get_filtration()) == list(st.get_filtration())
+
+    assert st.__is_persistence_defined() == True
+    assert st_copy.__is_persistence_defined() == False
+
+    st_copy.assign_filtration([1, 2, 3], 2.)
+    assert list(st_copy.get_filtration()) == list(st.get_filtration())
+
+    # test double free
+    del st
+    del st_copy
+
+def test_simplex_tree_deep_copy_constructor():
+    st = SimplexTree()
+    st.insert([1, 2, 3], 0.)
+    # persistence is not copied
+    st.compute_persistence()
+
+    st_copy = SimplexTree(st, copy = True)
+    # TODO(VR): when #463 is merged, replace with
+    # assert st_copy == st
+    assert st_copy.num_vertices() == st.num_vertices()
+    assert st_copy.num_simplices() == st.num_simplices()
+    st_filt_list = list(st.get_filtration())
+    assert list(st_copy.get_filtration()) == st_filt_list
+
+    assert st.__is_persistence_defined() == True
+    assert st_copy.__is_persistence_defined() == False
+
+    st_copy.remove_maximal_simplex([1, 2, 3])
+    a_filt_list = list(st_copy.get_filtration())
+    assert len(a_filt_list) < len(st_filt_list)
+
+    for a_splx in a_filt_list:
+        assert a_splx in st_filt_list
+    
+    # test double free
+    del st
+    del st_copy
+
+def test_simplex_tree_shallow_copy():
+    st = SimplexTree()
+    st.insert([1, 2, 3], 0.)
+    # persistence is not copied
+    st.compute_persistence()
+
+    st_copy = SimplexTree(st, copy = False)
+    # TODO(VR): when #463 is merged, replace with
+    # assert st_copy == st
+    assert st_copy.num_vertices() == st.num_vertices()
+    assert st_copy.num_simplices() == st.num_simplices()
+    assert list(st_copy.get_filtration()) == list(st.get_filtration())
+
+    assert st.__is_persistence_defined() == True
+    assert st_copy.__is_persistence_defined() == False
+
+    st_copy.assign_filtration([1, 2, 3], 2.)
+    assert list(st_copy.get_filtration()) == list(st.get_filtration())
+
+    # test double free
+    del st
+    del st_copy
-- 
cgit v1.2.3


From 43981a4d487669fe2002337ab62b72dd9e83a64a Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 14 Feb 2022 11:08:00 +0100
Subject: Remove shallow copy

---
 .github/next_release.md              |  2 +-
 src/python/gudhi/simplex_tree.pyx    | 54 +++++++++++++-----------------------
 src/python/test/test_simplex_tree.py | 50 ++-------------------------------
 3 files changed, 23 insertions(+), 83 deletions(-)

diff --git a/.github/next_release.md b/.github/next_release.md
index 3946404b..3d4761eb 100644
--- a/.github/next_release.md
+++ b/.github/next_release.md
@@ -14,7 +14,7 @@ Below is a list of changes made since GUDHI 3.5.0:
      - A more flexible Betti curve class capable of computing exact curves
 
 - [Simplex tree](https://gudhi.inria.fr/python/latest/simplex_tree_ref.html)
-     - `__copy__`, `__deepcopy__`, `copy` and copy constructors
+     - `__deepcopy__`, `copy` and copy constructors
 
 - Installation
      - Boost &ge; 1.66.0 is now required (was &ge; 1.56.0).
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index ed7c3b92..0213e363 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -30,7 +30,6 @@ cdef class SimplexTree:
     # unfortunately 'cdef public Simplex_tree_interface_full_featured* thisptr' is not possible
     # Use intptr_t instead to cast the pointer
     cdef public intptr_t thisptr
-    cdef bool __thisptr_to_be_deleted
 
     # Get the pointer casted as it should be
     cdef Simplex_tree_interface_full_featured* get_ptr(self) nogil:
@@ -39,36 +38,32 @@ cdef class SimplexTree:
     cdef Simplex_tree_persistence_interface * pcohptr
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, other = None, copy = True):
+    def __init__(self, other = None):
         """SimplexTree constructor.
-        :param other: If `other` is a SimplexTree (default = None), the SimplexTree is constructed from a deep/shallow copy of `other`.
+
+        :param other: If `other` is a `None` (default value), an empty `SimplexTree` is created.
+            If `other` is a `SimplexTree`, the `SimplexTree` is constructed from a deep copy of `other`.
         :type other: SimplexTree
-        :param copy: If `True`, the copy will be deep and if `False, the copy will be shallow. Default is `True`.
-        :type copy: bool
-        :returns: A simplex tree that is a (deep or shallow) copy of itself.
+        :returns: An empty or a copy simplex tree.
         :rtype: SimplexTree
-        :note: copy constructor requires :func:`compute_persistence` to be launched again as the result is not copied.
+
+        :note: If the `SimplexTree` is a copy, it requires :func:`compute_persistence` to be launched again as the
+            persistence result is not copied.
         """
 
     # The real cython constructor
-    def __cinit__(self, other = None, copy = True):
+    def __cinit__(self, other = None):
         cdef SimplexTree ostr
         if other and type(other) is SimplexTree:
             ostr = <SimplexTree> other
-            if copy:
-                self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured(dereference(ostr.get_ptr())))
-            else:
-                self.thisptr = ostr.thisptr
-                # Avoid double free - The original is in charge of deletion
-                self.__thisptr_to_be_deleted = False
+            self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured(dereference(ostr.get_ptr())))
         else:
-            self.__thisptr_to_be_deleted = True
             self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured())
 
     def __dealloc__(self):
         cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr()
         # Avoid double free - The original is in charge of deletion
-        if ptr != NULL and self.__thisptr_to_be_deleted:
+        if ptr != NULL:
             del ptr
         if self.pcohptr != NULL:
             del self.pcohptr
@@ -83,33 +78,24 @@ cdef class SimplexTree:
          """
         return self.pcohptr != NULL
 
-    def copy(self, deep=True):
+    def copy(self):
         """ 
-        :param deep: If `True`, the copy will be deep and if `False`, the copy will be shallow. Default is `True`.
-        :type deep: bool
-        :returns: A simplex tree that is a (deep or shallow) copy of itself.
+        :returns: A simplex tree that is a deep copy of itself.
         :rtype: SimplexTree
-        :note: copy requires :func:`compute_persistence` to be launched again as the result is not copied.
+
+        :note: copy requires :func:`compute_persistence` to be launched again as the persistence result is not copied.
         """
         stree = SimplexTree()
         cdef Simplex_tree_interface_full_featured* stree_ptr
         cdef Simplex_tree_interface_full_featured* self_ptr=self.get_ptr()
-        if deep:
-            with nogil:
-                stree_ptr = new Simplex_tree_interface_full_featured(dereference(self_ptr))
-    
-            stree.thisptr = <intptr_t>(stree_ptr)
-        else:
-            stree.thisptr = self.thisptr
-            # Avoid double free - The original is in charge of deletion
-            stree.__thisptr_to_be_deleted = False
-        return stree
+        with nogil:
+            stree_ptr = new Simplex_tree_interface_full_featured(dereference(self_ptr))
 
-    def __copy__(self):
-        return self.copy(deep=False)
+        stree.thisptr = <intptr_t>(stree_ptr)
+        return stree
 
     def __deepcopy__(self):
-        return self.copy(deep=True)
+        return self.copy()
 
     def filtration(self, simplex):
         """This function returns the filtration value for a given N-simplex in
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 6db6d8fb..62dcc865 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -454,7 +454,7 @@ def test_simplex_tree_deep_copy():
     # persistence is not copied
     st.compute_persistence()
 
-    st_copy = st.copy(deep=True)
+    st_copy = st.copy()
     # TODO(VR): when #463 is merged, replace with
     # assert st_copy == st
     assert st_copy.num_vertices() == st.num_vertices()
@@ -476,36 +476,13 @@ def test_simplex_tree_deep_copy():
     del st
     del st_copy
 
-def test_simplex_tree_shallow_copy():
-    st = SimplexTree()
-    st.insert([1, 2, 3], 0.)
-    # persistence is not copied
-    st.compute_persistence()
-
-    st_copy = st.copy(deep=False)
-    # TODO(VR): when #463 is merged, replace with
-    # assert st_copy == st
-    assert st_copy.num_vertices() == st.num_vertices()
-    assert st_copy.num_simplices() == st.num_simplices()
-    assert list(st_copy.get_filtration()) == list(st.get_filtration())
-
-    assert st.__is_persistence_defined() == True
-    assert st_copy.__is_persistence_defined() == False
-
-    st_copy.assign_filtration([1, 2, 3], 2.)
-    assert list(st_copy.get_filtration()) == list(st.get_filtration())
-
-    # test double free
-    del st
-    del st_copy
-
 def test_simplex_tree_deep_copy_constructor():
     st = SimplexTree()
     st.insert([1, 2, 3], 0.)
     # persistence is not copied
     st.compute_persistence()
 
-    st_copy = SimplexTree(st, copy = True)
+    st_copy = SimplexTree(st)
     # TODO(VR): when #463 is merged, replace with
     # assert st_copy == st
     assert st_copy.num_vertices() == st.num_vertices()
@@ -526,26 +503,3 @@ def test_simplex_tree_deep_copy_constructor():
     # test double free
     del st
     del st_copy
-
-def test_simplex_tree_shallow_copy():
-    st = SimplexTree()
-    st.insert([1, 2, 3], 0.)
-    # persistence is not copied
-    st.compute_persistence()
-
-    st_copy = SimplexTree(st, copy = False)
-    # TODO(VR): when #463 is merged, replace with
-    # assert st_copy == st
-    assert st_copy.num_vertices() == st.num_vertices()
-    assert st_copy.num_simplices() == st.num_simplices()
-    assert list(st_copy.get_filtration()) == list(st.get_filtration())
-
-    assert st.__is_persistence_defined() == True
-    assert st_copy.__is_persistence_defined() == False
-
-    st_copy.assign_filtration([1, 2, 3], 2.)
-    assert list(st_copy.get_filtration()) == list(st.get_filtration())
-
-    # test double free
-    del st
-    del st_copy
-- 
cgit v1.2.3


From c7303733b28d3bc429d4bdfd030409e07430599d Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 14 Feb 2022 11:29:49 +0100
Subject: Rewrite test as #463 is merged. Rephrase constructor documentation

---
 src/python/gudhi/simplex_tree.pyx    |  4 ++--
 src/python/test/test_simplex_tree.py | 20 ++++++++------------
 2 files changed, 10 insertions(+), 14 deletions(-)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 6701d98d..a2b8719a 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -41,9 +41,9 @@ cdef class SimplexTree:
     def __init__(self, other = None):
         """SimplexTree constructor.
 
-        :param other: If `other` is a `None` (default value), an empty `SimplexTree` is created.
+        :param other: If `other` is `None` (default value), an empty `SimplexTree` is created.
             If `other` is a `SimplexTree`, the `SimplexTree` is constructed from a deep copy of `other`.
-        :type other: SimplexTree
+        :type other: SimplexTree (Optional)
         :returns: An empty or a copy simplex tree.
         :rtype: SimplexTree
 
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 32fc63ec..2c2e09a2 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -463,20 +463,18 @@ def test_equality_operator():
 def test_simplex_tree_deep_copy():
     st = SimplexTree()
     st.insert([1, 2, 3], 0.)
-    # persistence is not copied
+    # compute persistence only on the original
     st.compute_persistence()
 
     st_copy = st.copy()
-    # TODO(VR): when #463 is merged, replace with
-    # assert st_copy == st
-    assert st_copy.num_vertices() == st.num_vertices()
-    assert st_copy.num_simplices() == st.num_simplices()
+    assert st_copy == st
     st_filt_list = list(st.get_filtration())
-    assert list(st_copy.get_filtration()) == st_filt_list
 
+    # check persistence is not copied
     assert st.__is_persistence_defined() == True
     assert st_copy.__is_persistence_defined() == False
 
+    # remove something in the copy and check the copy is included in the original
     st_copy.remove_maximal_simplex([1, 2, 3])
     a_filt_list = list(st_copy.get_filtration())
     assert len(a_filt_list) < len(st_filt_list)
@@ -491,20 +489,18 @@ def test_simplex_tree_deep_copy():
 def test_simplex_tree_deep_copy_constructor():
     st = SimplexTree()
     st.insert([1, 2, 3], 0.)
-    # persistence is not copied
+    # compute persistence only on the original
     st.compute_persistence()
 
     st_copy = SimplexTree(st)
-    # TODO(VR): when #463 is merged, replace with
-    # assert st_copy == st
-    assert st_copy.num_vertices() == st.num_vertices()
-    assert st_copy.num_simplices() == st.num_simplices()
+    assert st_copy == st
     st_filt_list = list(st.get_filtration())
-    assert list(st_copy.get_filtration()) == st_filt_list
 
+    # check persistence is not copied
     assert st.__is_persistence_defined() == True
     assert st_copy.__is_persistence_defined() == False
 
+    # remove something in the copy and check the copy is included in the original
     st_copy.remove_maximal_simplex([1, 2, 3])
     a_filt_list = list(st_copy.get_filtration())
     assert len(a_filt_list) < len(st_filt_list)
-- 
cgit v1.2.3


From 1a0551c033e55024cd0f00302cd9df1f356bab44 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 14 Feb 2022 11:31:53 +0100
Subject: some left over

---
 src/python/gudhi/simplex_tree.pyx | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index a2b8719a..8f760422 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -62,7 +62,6 @@ cdef class SimplexTree:
 
     def __dealloc__(self):
         cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr()
-        # Avoid double free - The original is in charge of deletion
         if ptr != NULL:
             del ptr
         if self.pcohptr != NULL:
-- 
cgit v1.2.3


From d5ac245a6dc4ab2d6e30689fc5d95503c40b6187 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 14 Feb 2022 12:31:30 +0100
Subject: code review: ctor shall raise a TypeError when constructed from
 something else than a SimplexTree

---
 src/python/gudhi/simplex_tree.pyx    | 10 +++++++---
 src/python/test/test_simplex_tree.py |  4 ++++
 2 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 8f760422..e1685ded 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -47,6 +47,7 @@ cdef class SimplexTree:
         :returns: An empty or a copy simplex tree.
         :rtype: SimplexTree
 
+        :raises TypeError: In case `other` is neither `None`, nor a `SimplexTree`.
         :note: If the `SimplexTree` is a copy, it requires :func:`compute_persistence` to be launched again as the
             persistence result is not copied.
         """
@@ -54,9 +55,12 @@ cdef class SimplexTree:
     # The real cython constructor
     def __cinit__(self, other = None):
         cdef SimplexTree ostr
-        if other and type(other) is SimplexTree:
-            ostr = <SimplexTree> other
-            self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured(dereference(ostr.get_ptr())))
+        if other:
+            if type(other) is SimplexTree:
+                ostr = <SimplexTree> other
+                self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured(dereference(ostr.get_ptr())))
+            else:
+                raise TypeError("`other` argument requires to be of type `SimplexTree`, or `None`.")
         else:
             self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured())
 
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 2c2e09a2..a8180ce8 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -511,3 +511,7 @@ def test_simplex_tree_deep_copy_constructor():
     # test double free
     del st
     del st_copy
+
+def test_simplex_tree_constructor_exception():
+    with pytest.raises(TypeError):
+        st = SimplexTree(other = "Construction from a string shall raise an exception")
-- 
cgit v1.2.3


From 4e63e52758daa7ef27782b4f129b7e55fa73de43 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 28 Feb 2022 09:50:08 +0100
Subject: code review: use 'isinstance' instead of 'type'

---
 src/python/gudhi/simplex_tree.pyx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index e1685ded..711796d4 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -56,7 +56,7 @@ cdef class SimplexTree:
     def __cinit__(self, other = None):
         cdef SimplexTree ostr
         if other:
-            if type(other) is SimplexTree:
+            if isinstance(other, SimplexTree):
                 ostr = <SimplexTree> other
                 self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured(dereference(ostr.get_ptr())))
             else:
-- 
cgit v1.2.3


From e7ba8967e9b5c3c7260522003d8f87e643b7912e Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 28 Feb 2022 10:56:47 +0100
Subject: doc review: reformulate the copy notes

---
 src/python/gudhi/simplex_tree.pyx | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 711796d4..91e079aa 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -48,8 +48,8 @@ cdef class SimplexTree:
         :rtype: SimplexTree
 
         :raises TypeError: In case `other` is neither `None`, nor a `SimplexTree`.
-        :note: If the `SimplexTree` is a copy, it requires :func:`compute_persistence` to be launched again as the
-            persistence result is not copied.
+        :note: If the `SimplexTree` is a copy, the persistence information is not copied. If you need it in the clone,
+            you have to call :func:`compute_persistence` on it even if you had already computed it in the original.
         """
 
     # The real cython constructor
@@ -86,7 +86,8 @@ cdef class SimplexTree:
         :returns: A simplex tree that is a deep copy of itself.
         :rtype: SimplexTree
 
-        :note: copy requires :func:`compute_persistence` to be launched again as the persistence result is not copied.
+        :note: The persistence information is not copied. If you need it in the clone, you have to call
+            :func:`compute_persistence` on it even if you had already computed it in the original.
         """
         stree = SimplexTree()
         cdef Simplex_tree_interface_full_featured* stree_ptr
-- 
cgit v1.2.3


From 41a976cc85ab36dc2df748b9d6900d77e76b2fa7 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 28 Feb 2022 14:58:24 +0100
Subject: code review: copy simplex tree c++ pointer in a factorized nogil
 function

---
 src/python/gudhi/simplex_tree.pyx | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 91e079aa..b8fabf78 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -54,11 +54,9 @@ cdef class SimplexTree:
 
     # The real cython constructor
     def __cinit__(self, other = None):
-        cdef SimplexTree ostr
         if other:
             if isinstance(other, SimplexTree):
-                ostr = <SimplexTree> other
-                self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured(dereference(ostr.get_ptr())))
+                self.thisptr = _get_copy_intptr(other)
             else:
                 raise TypeError("`other` argument requires to be of type `SimplexTree`, or `None`.")
         else:
@@ -90,12 +88,7 @@ cdef class SimplexTree:
             :func:`compute_persistence` on it even if you had already computed it in the original.
         """
         stree = SimplexTree()
-        cdef Simplex_tree_interface_full_featured* stree_ptr
-        cdef Simplex_tree_interface_full_featured* self_ptr=self.get_ptr()
-        with nogil:
-            stree_ptr = new Simplex_tree_interface_full_featured(dereference(self_ptr))
-
-        stree.thisptr = <intptr_t>(stree_ptr)
+        stree.thisptr = _get_copy_intptr(self)
         return stree
 
     def __deepcopy__(self):
@@ -687,3 +680,6 @@ cdef class SimplexTree:
         :rtype: bool
         """
         return dereference(self.get_ptr()) == dereference(other.get_ptr())
+
+cdef intptr_t _get_copy_intptr(SimplexTree stree) nogil:
+    return <intptr_t>(new Simplex_tree_interface_full_featured(dereference(stree.get_ptr())))
-- 
cgit v1.2.3


From 533612bcd73c2aa8c1e9e2209973a0e999ab3f41 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 7 Mar 2022 09:21:36 +0100
Subject: tangential - include simplex tree when used

---
 src/Tangential_complex/benchmark/benchmark_tc.cpp       | 1 +
 src/Tangential_complex/example/example_basic.cpp        | 4 +---
 src/Tangential_complex/example/example_with_perturb.cpp | 1 +
 src/Tangential_complex/test/test_tangential_complex.cpp | 1 +
 4 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/Tangential_complex/benchmark/benchmark_tc.cpp b/src/Tangential_complex/benchmark/benchmark_tc.cpp
index 6da1425f..8e7c72ff 100644
--- a/src/Tangential_complex/benchmark/benchmark_tc.cpp
+++ b/src/Tangential_complex/benchmark/benchmark_tc.cpp
@@ -33,6 +33,7 @@ const std::size_t ONLY_LOAD_THE_FIRST_N_POINTS = 20000000;
 #include <gudhi/sparsify_point_set.h>
 #include <gudhi/random_point_generators.h>
 #include <gudhi/Tangential_complex/utilities.h>
+#include <gudhi/Simplex_tree.h>
 
 #include <CGAL/assertions_behaviour.h>
 #include <CGAL/Epick_d.h>
diff --git a/src/Tangential_complex/example/example_basic.cpp b/src/Tangential_complex/example/example_basic.cpp
index ab35edf0..c50b9b8c 100644
--- a/src/Tangential_complex/example/example_basic.cpp
+++ b/src/Tangential_complex/example/example_basic.cpp
@@ -1,7 +1,6 @@
 #include <gudhi/Tangential_complex.h>
 #include <gudhi/sparsify_point_set.h>
-//#include <gudhi/Fake_simplex_tree.h>
-
+#include <gudhi/Simplex_tree.h>
 
 #include <CGAL/Epick_d.h>
 #include <CGAL/Random.h>
@@ -39,7 +38,6 @@ int main(void) {
 
   // Export the TC into a Simplex_tree
   Gudhi::Simplex_tree<> stree;
-  //Gudhi::Fake_simplex_tree stree;
   tc.create_complex(stree);
 
   // Display stats about inconsistencies
diff --git a/src/Tangential_complex/example/example_with_perturb.cpp b/src/Tangential_complex/example/example_with_perturb.cpp
index d0d877ea..e70e2980 100644
--- a/src/Tangential_complex/example/example_with_perturb.cpp
+++ b/src/Tangential_complex/example/example_with_perturb.cpp
@@ -1,5 +1,6 @@
 #include <gudhi/Tangential_complex.h>
 #include <gudhi/sparsify_point_set.h>
+#include <gudhi/Simplex_tree.h>
 
 #include <CGAL/Epick_d.h>
 #include <CGAL/Random.h>
diff --git a/src/Tangential_complex/test/test_tangential_complex.cpp b/src/Tangential_complex/test/test_tangential_complex.cpp
index 023c1e1a..a24b9ae2 100644
--- a/src/Tangential_complex/test/test_tangential_complex.cpp
+++ b/src/Tangential_complex/test/test_tangential_complex.cpp
@@ -14,6 +14,7 @@
 
 #include <gudhi/Tangential_complex.h>
 #include <gudhi/sparsify_point_set.h>
+#include <gudhi/Simplex_tree.h>
 
 #include <CGAL/Epick_d.h>
 #include <CGAL/Random.h>
-- 
cgit v1.2.3


From d114eaf05a4d622045acea1d7fdced20e7139e53 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 21 Mar 2022 16:09:05 +0100
Subject: code review: assign_children(new_sib) was happening too late

---
 src/Simplex_tree/include/gudhi/Simplex_tree.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 85790baf..d48f6616 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -1283,6 +1283,7 @@ class Simplex_tree {
         Siblings * new_sib = new Siblings(siblings,  // oncles
                                           simplex->first,  // parent
                                           boost::adaptors::reverse(intersection));  // boost::container::ordered_unique_range_t
+        simplex->second.assign_children(new_sib);
         std::vector<Vertex_handle> blocked_new_sib_vertex_list;
         // As all intersections are inserted, we can call the blocker function on all new_sib members
         for (auto new_sib_member = new_sib->members().begin();
@@ -1305,7 +1306,6 @@ class Simplex_tree {
             new_sib->members().erase(blocked_new_sib_member);
           }
           // ensure recursive call
-          simplex->second.assign_children(new_sib);
           siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex);
         }
       } else {
-- 
cgit v1.2.3


From 8667b75601cc0420f3b196d5b44ec93fd6058057 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 21 Mar 2022 17:14:30 +0100
Subject: Add some c++ tests and use unitary_tests_utils. Fix python test

---
 .../simplex_tree_graph_expansion_unit_test.cpp     | 192 ++++++++++++++++-----
 src/python/test/test_simplex_tree.py               |  11 +-
 2 files changed, 151 insertions(+), 52 deletions(-)

diff --git a/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp
index 881a06ae..6d63d8ae 100644
--- a/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp
+++ b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp
@@ -9,33 +9,62 @@
  */
 
 #include <iostream>
-#include <fstream>
-#include <string>
-#include <algorithm>
-#include <utility> // std::pair, std::make_pair
-#include <cmath> // float comparison
-#include <limits>
-#include <functional> // greater
+#include <vector>
 
 #define BOOST_TEST_DYN_LINK
-#define BOOST_TEST_MODULE "simplex_tree"
+#define BOOST_TEST_MODULE "simplex_tree_graph_expansion"
 #include <boost/test/unit_test.hpp>
 #include <boost/mpl/list.hpp>
 
-//  ^
-// /!\ Nothing else from Simplex_tree shall be included to test includes are well defined.
 #include "gudhi/Simplex_tree.h"
+#include <gudhi/Unitary_tests_utils.h>
 
 using namespace Gudhi;
 
 typedef boost::mpl::list<Simplex_tree<>, Simplex_tree<Simplex_tree_options_fast_persistence>> list_of_tested_variants;
 
+BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_all_is_blocked, typeST, list_of_tested_variants) {
+  std::clog << "********************************************************************\n";
+  std::clog << "simplex_tree_expansion_all_is_blocked\n";
+  std::clog << "********************************************************************\n";
+  using Simplex_handle = typename typeST::Simplex_handle;
+  // Construct the Simplex Tree with a 1-skeleton graph example
+  typeST simplex_tree;
+
+  simplex_tree.insert_simplex({0, 1}, 0.);
+  simplex_tree.insert_simplex({0, 2}, 1.);
+  simplex_tree.insert_simplex({0, 3}, 2.);
+  simplex_tree.insert_simplex({1, 2}, 3.);
+  simplex_tree.insert_simplex({1, 3}, 4.);
+  simplex_tree.insert_simplex({2, 3}, 5.);
+  simplex_tree.insert_simplex({2, 4}, 6.);
+  simplex_tree.insert_simplex({3, 6}, 7.);
+  simplex_tree.insert_simplex({4, 5}, 8.);
+  simplex_tree.insert_simplex({4, 6}, 9.);
+  simplex_tree.insert_simplex({5, 6}, 10.);
+  simplex_tree.insert_simplex({6}, 10.);
 
-bool AreAlmostTheSame(float a, float b) {
-  return std::fabs(a - b) < std::numeric_limits<float>::epsilon();
+  typeST stree_copy = simplex_tree;
+
+  simplex_tree.expansion_with_blockers(3, [&](Simplex_handle sh){ return true; });
+
+  std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices";
+  std::clog << " - dimension " << simplex_tree.dimension() << "\n";
+  std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
+  for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
+    std::clog << "   " << "[" << simplex_tree.filtration(f_simplex) << "] ";
+    for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex))
+      std::clog << "(" << vertex << ")";
+    std::clog << std::endl;
+  }
+
+  BOOST_CHECK(stree_copy == simplex_tree);
 }
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_3, typeST, list_of_tested_variants) {
+  std::clog << "********************************************************************\n";
+  std::clog << "simplex_tree_expansion_with_blockers_3\n";
+  std::clog << "********************************************************************\n";
   using Simplex_handle = typename typeST::Simplex_handle;
   // Construct the Simplex Tree with a 1-skeleton graph example
   typeST simplex_tree;
@@ -72,9 +101,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_3, typeST, li
       return result;
     });
 
-  std::clog << "********************************************************************\n";
-  std::clog << "simplex_tree_expansion_with_blockers_3\n";
-  std::clog << "********************************************************************\n";
   std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices";
   std::clog << " - dimension " << simplex_tree.dimension() << "\n";
   std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
@@ -89,15 +115,23 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_3, typeST, li
   BOOST_CHECK(simplex_tree.dimension() == 3);
   // {4, 5, 6} shall be blocked
   BOOST_CHECK(simplex_tree.find({4, 5, 6}) == simplex_tree.null_simplex());
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2})), 4.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,3})), 5.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,2,3})), 6.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({1,2,3})), 6.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2,3})), 7.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2})),
+                                                          static_cast<typename typeST::Filtration_value>(4.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,3})),
+                                                          static_cast<typename typeST::Filtration_value>(5.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(6.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({1,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(6.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(7.));
 
 }
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_2, typeST, list_of_tested_variants) {
+  std::clog << "********************************************************************\n";
+  std::clog << "simplex_tree_expansion_with_blockers_2\n";
+  std::clog << "********************************************************************\n";
   using Simplex_handle = typename typeST::Simplex_handle;
   // Construct the Simplex Tree with a 1-skeleton graph example
   typeST simplex_tree;
@@ -134,9 +168,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_2, typeST, li
       return result;
     });
 
-  std::clog << "********************************************************************\n";
-  std::clog << "simplex_tree_expansion_with_blockers_2\n";
-  std::clog << "********************************************************************\n";
   std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices";
   std::clog << " - dimension " << simplex_tree.dimension() << "\n";
   std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
@@ -151,14 +182,22 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_2, typeST, li
   BOOST_CHECK(simplex_tree.dimension() == 2);
   // {4, 5, 6} shall be blocked
   BOOST_CHECK(simplex_tree.find({4, 5, 6}) == simplex_tree.null_simplex());
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2})), 4.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,3})), 5.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,2,3})), 6.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({1,2,3})), 6.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2})),
+                                                          static_cast<typename typeST::Filtration_value>(4.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,3})),
+                                                          static_cast<typename typeST::Filtration_value>(5.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(6.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({1,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(6.));
   BOOST_CHECK(simplex_tree.find({0,1,2,3}) == simplex_tree.null_simplex());
 }
 
-BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion, typeST, list_of_tested_variants) {
+BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_find_simplex_blockers, typeST, list_of_tested_variants) {
+  std::clog << "********************************************************************\n";
+  std::clog << "simplex_tree_expansion_with_find_simplex_blockers\n";
+  std::clog << "********************************************************************\n";
+  using Simplex_handle = typename typeST::Simplex_handle;
   // Construct the Simplex Tree with a 1-skeleton graph example
   typeST simplex_tree;
 
@@ -175,10 +214,66 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion, typeST, list_of_tested_var
   simplex_tree.insert_simplex({5, 6}, 10.);
   simplex_tree.insert_simplex({6}, 10.);
 
-  simplex_tree.expansion(3);
+  simplex_tree.expansion_with_blockers(3, [&](Simplex_handle sh){
+      bool result = false;
+      std::clog << "Blocker on [";
+      std::vector<typename typeST::Vertex_handle> simplex;
+      // User can loop on the vertices from the given simplex_handle i.e.
+      for (auto vertex : simplex_tree.simplex_vertex_range(sh)) {
+        // We block the expansion, if the vertex '1' is in the given list of vertices
+        if (vertex == 1)
+          result = true;
+        std::clog << vertex << ", ";
+        simplex.push_back(vertex);
+      }
+      std::clog << "] => " << result << std::endl;
+      // Not efficient but test it works - required by the python interface
+      BOOST_CHECK(simplex_tree.find(simplex) == sh);
+      return result;
+    });
+
+  std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices";
+  std::clog << " - dimension " << simplex_tree.dimension() << "\n";
+  std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
+  for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
+    std::clog << "   " << "[" << simplex_tree.filtration(f_simplex) << "] ";
+    for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex))
+      std::clog << "(" << vertex << ")";
+    std::clog << std::endl;
+  }
+
+  BOOST_CHECK(simplex_tree.num_simplices() == 20);
+  BOOST_CHECK(simplex_tree.dimension() == 2);
+
+  // {1, 2, 3}, {0, 1, 2} and {0, 1, 3} shall be blocked as it contains vertex 1
+  BOOST_CHECK(simplex_tree.find({4, 5, 6}) != simplex_tree.null_simplex());
+  BOOST_CHECK(simplex_tree.find({1, 2, 3}) == simplex_tree.null_simplex());
+  BOOST_CHECK(simplex_tree.find({0, 2, 3}) != simplex_tree.null_simplex());
+  BOOST_CHECK(simplex_tree.find({0, 1, 2}) == simplex_tree.null_simplex());
+  BOOST_CHECK(simplex_tree.find({0, 1, 3}) == simplex_tree.null_simplex());
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_3, typeST, list_of_tested_variants) {
   std::clog << "********************************************************************\n";
   std::clog << "simplex_tree_expansion_3\n";
   std::clog << "********************************************************************\n";
+  // Construct the Simplex Tree with a 1-skeleton graph example
+  typeST simplex_tree;
+
+  simplex_tree.insert_simplex({0, 1}, 0.);
+  simplex_tree.insert_simplex({0, 2}, 1.);
+  simplex_tree.insert_simplex({0, 3}, 2.);
+  simplex_tree.insert_simplex({1, 2}, 3.);
+  simplex_tree.insert_simplex({1, 3}, 4.);
+  simplex_tree.insert_simplex({2, 3}, 5.);
+  simplex_tree.insert_simplex({2, 4}, 6.);
+  simplex_tree.insert_simplex({3, 6}, 7.);
+  simplex_tree.insert_simplex({4, 5}, 8.);
+  simplex_tree.insert_simplex({4, 6}, 9.);
+  simplex_tree.insert_simplex({5, 6}, 10.);
+  simplex_tree.insert_simplex({6}, 10.);
+
+  simplex_tree.expansion(3);
   std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices";
   std::clog << " - dimension " << simplex_tree.dimension() << "\n";
   std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
@@ -192,16 +287,25 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion, typeST, list_of_tested_var
   BOOST_CHECK(simplex_tree.num_simplices() == 24);
   BOOST_CHECK(simplex_tree.dimension() == 3);
 
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({4,5,6})), 10.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2})), 3.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,3})), 4.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,2,3})), 5.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({1,2,3})), 5.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2,3})), 5.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({4,5,6})),
+                                                          static_cast<typename typeST::Filtration_value>(10.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2})),
+                                                          static_cast<typename typeST::Filtration_value>(3.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,3})),
+                                                          static_cast<typename typeST::Filtration_value>(4.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(5.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({1,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(5.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(5.));
 
 }
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_2, typeST, list_of_tested_variants) {
+  std::clog << "********************************************************************\n";
+  std::clog << "simplex_tree_expansion_2\n";
+  std::clog << "********************************************************************\n";
   // Construct the Simplex Tree with a 1-skeleton graph example
   typeST simplex_tree;
 
@@ -220,9 +324,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_2, typeST, list_of_tested_v
 
   simplex_tree.expansion(2);
 
-  std::clog << "********************************************************************\n";
-  std::clog << "simplex_tree_expansion_2\n";
-  std::clog << "********************************************************************\n";
   std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices";
   std::clog << " - dimension " << simplex_tree.dimension() << "\n";
   std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
@@ -236,10 +337,15 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_2, typeST, list_of_tested_v
   BOOST_CHECK(simplex_tree.num_simplices() == 23);
   BOOST_CHECK(simplex_tree.dimension() == 2);
 
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({4,5,6})), 10.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2})), 3.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,3})), 4.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,2,3})), 5.));
-  BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({1,2,3})), 5.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({4,5,6})),
+                                                          static_cast<typename typeST::Filtration_value>(10.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2})),
+                                                          static_cast<typename typeST::Filtration_value>(3.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,3})),
+                                                          static_cast<typename typeST::Filtration_value>(4.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(5.));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({1,2,3})),
+                                                          static_cast<typename typeST::Filtration_value>(5.));
   BOOST_CHECK(simplex_tree.find({0,1,2,3}) == simplex_tree.null_simplex());
 }
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index d8173b52..eb481a49 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -531,8 +531,6 @@ def test_expansion_with_blocker():
     st.insert([5,6],10)
     st.insert([6],10)
 
-    # One cannot modify filtration inside blocker - segfault - Let's record accepted simplices
-    accepted_simp = []
     def blocker(simplex):
         try:
             # Block all simplices that countains vertex 6
@@ -541,12 +539,10 @@ def test_expansion_with_blocker():
             return True
         except ValueError:
             print(simplex, ' is accepted')
-            accepted_simp.append(simplex)
+            st.assign_filtration(simplex, st.filtration(simplex) + 1.)
             return False
 
     st.expansion_with_blocker(2, blocker)
-    for simplex in accepted_simp:
-        st.assign_filtration(simplex, st.filtration(simplex) + 1.)
     assert st.num_simplices() == 22
     assert st.dimension() == 2
     assert st.find([4,5,6]) == False
@@ -555,10 +551,7 @@ def test_expansion_with_blocker():
     assert st.filtration([0,2,3]) == 6.
     assert st.filtration([1,2,3]) == 6.
 
-    accepted_simp = []
     st.expansion_with_blocker(3, blocker)
-    for simplex in accepted_simp:
-        st.assign_filtration(simplex, st.filtration(simplex) + 1.)
     assert st.num_simplices() == 23
     assert st.dimension() == 3
     assert st.find([4,5,6]) == False
@@ -566,4 +559,4 @@ def test_expansion_with_blocker():
     assert st.filtration([0,1,3]) == 5.
     assert st.filtration([0,2,3]) == 6.
     assert st.filtration([1,2,3]) == 6.
-    assert st.filtration([0,1,2,3]) == 6.
+    assert st.filtration([0,1,2,3]) == 7.
-- 
cgit v1.2.3


From 924302f092d6d08c59500c5f6e75bcad9416581a Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 28 Mar 2022 08:39:21 +0200
Subject: doc review: remove segfault note as it is no more the case

---
 src/python/gudhi/simplex_tree.pyx | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 644110c1..fc704358 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -487,9 +487,6 @@ cdef class SimplexTree:
         with this new simplex (represented as a list of int). If `blocker_func` returns `True`, the simplex is removed,
         otherwise it is kept. The algorithm then proceeds with the next candidate.
 
-        Note that you cannot update the filtration value of the simplex during the evaluation of `blocker_func`, as it
-        would segfault. 
-
         :param max_dim: Expansion maximal dimension value.
         :type max_dim: int
         :param blocker_func: Blocker oracle.
-- 
cgit v1.2.3


From f6a45247e9a9f126c214d1b1003ae19fb2cc84a3 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 28 Mar 2022 09:00:23 +0200
Subject: doc review: add a warning about phantom simplices for blocker_func

---
 src/python/gudhi/simplex_tree.pyx | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index fc704358..a4914184 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -487,6 +487,11 @@ cdef class SimplexTree:
         with this new simplex (represented as a list of int). If `blocker_func` returns `True`, the simplex is removed,
         otherwise it is kept. The algorithm then proceeds with the next candidate.
 
+        .. warning::
+            Several candidates of the same dimension may be inserted simultaneously before calling `block_simplex`, so
+            if you examine the complex in `block_simplex`, you may hit a few simplices of the same dimension that have
+            not been vetted by `block_simplex` yet, or have already been rejected but not yet removed.
+
         :param max_dim: Expansion maximal dimension value.
         :type max_dim: int
         :param blocker_func: Blocker oracle.
-- 
cgit v1.2.3


From 53818204ba2f914c89e21deb6bfe7cb4a0122718 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Wed, 30 Mar 2022 09:22:09 +0200
Subject: Remove this useless line

---
 src/Alpha_complex/include/gudhi/Alpha_complex_3d.h | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
index df5c630e..b3dbc9bb 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
@@ -57,8 +57,6 @@ namespace Gudhi {
 
 namespace alpha_complex {
 
-thread_local double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001;
-
 // Value_from_iterator returns the filtration value from an iterator on alpha shapes values
 //
 // FAST                         SAFE                         EXACT
-- 
cgit v1.2.3


From 495f3d9b51f5ed48151118fedcc61f5067f04d51 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Wed, 30 Mar 2022 12:22:56 +0200
Subject: get/set float_relative_precision and its documentation

---
 src/python/doc/alpha_complex_user.rst        |  3 ++-
 src/python/gudhi/alpha_complex.pyx           | 28 ++++++++++++++++++++++++++++
 src/python/include/Alpha_complex_interface.h | 10 ++++++++++
 3 files changed, 40 insertions(+), 1 deletion(-)

diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst
index cfd22742..b060c86e 100644
--- a/src/python/doc/alpha_complex_user.rst
+++ b/src/python/doc/alpha_complex_user.rst
@@ -27,7 +27,8 @@ Remarks
   If you pass :code:`precision = 'exact'` to :func:`~gudhi.AlphaComplex.__init__`, the filtration values are the exact
   ones converted to float. This can be very slow.
   If you pass :code:`precision = 'safe'` (the default), the filtration values are only
-  guaranteed to have a small multiplicative error compared to the exact value.
+  guaranteed to have a small multiplicative error compared to the exact value, see
+  :func:`~gudhi.AlphaComplex.set_float_relative_precision` to modify the precision.
   A drawback, when computing persistence, is that an empty exact interval [10^12,10^12] may become a
   non-empty approximate interval [10^12,10^12+10^6].
   Using :code:`precision = 'fast'` makes the computations slightly faster, and the combinatorics are still exact, but
diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index a4888914..9af62b43 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -31,6 +31,10 @@ cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
         Alpha_complex_interface(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except +
         vector[double] get_point(int vertex) nogil except +
         void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except +
+        @staticmethod
+        void set_float_relative_precision(double precision) nogil
+        @staticmethod
+        double get_float_relative_precision() nogil
 
 # AlphaComplex python interface
 cdef class AlphaComplex:
@@ -133,3 +137,27 @@ cdef class AlphaComplex:
             self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr,
                                               mas, compute_filtration)
         return stree
+
+    @staticmethod
+    def set_float_relative_precision(precision):
+        """
+        :param precision: When the AlphaComplex is constructed with :code:`precision = 'safe'` (the default),
+            one can set the float relative precision of filtration values computed in
+            :func:`~gudhi.AlphaComplex.create_simplex_tree`.
+            Default is :code:`1e-5` (cf. :func:`~gudhi.AlphaComplex.get_float_relative_precision`).
+            For more details, please refer to
+            `CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double <https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html>`_
+        :type precision: float
+        """
+        if precision <=0. or precision >= 1.:
+            raise ValueError("Relative precision value must be strictly greater than 0 and strictly lower than 1")
+        Alpha_complex_interface.set_float_relative_precision(precision)
+    
+    @staticmethod
+    def get_float_relative_precision():
+        """
+        :returns: The float relative precision of filtration values computation in
+            :func:`~gudhi.AlphaComplex.create_simplex_tree`.
+        :rtype: float
+        """
+        return Alpha_complex_interface.get_float_relative_precision()
diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h
index 671af4a4..469b91ce 100644
--- a/src/python/include/Alpha_complex_interface.h
+++ b/src/python/include/Alpha_complex_interface.h
@@ -57,6 +57,16 @@ class Alpha_complex_interface {
       alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value);
   }
 
+  static void set_float_relative_precision(double precision) {
+    // cf. Exact_alpha_complex_dD kernel type in Alpha_complex_factory.h
+    CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>::FT::set_relative_precision_of_to_double(precision);
+  }
+
+  static double get_float_relative_precision() {
+    // cf. Exact_alpha_complex_dD kernel type in Alpha_complex_factory.h
+    return CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>::FT::get_relative_precision_of_to_double();
+  }
+
  private:
   std::unique_ptr<Abstract_alpha_complex> alpha_ptr_;
 };
-- 
cgit v1.2.3


From 37dcedb59cd2a6b70119688ea24c978749d63154 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Thu, 31 Mar 2022 16:37:20 +0200
Subject: Test the feature

---
 src/python/test/test_alpha_complex.py | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py
index f15284f3..f81e6137 100755
--- a/src/python/test/test_alpha_complex.py
+++ b/src/python/test/test_alpha_complex.py
@@ -286,3 +286,30 @@ def _weighted_doc_example(precision):
 def test_weighted_doc_example():
     for precision in ['fast', 'safe', 'exact']:
         _weighted_doc_example(precision)
+
+def test_float_relative_precision():
+    assert AlphaComplex.get_float_relative_precision() == 1e-5
+    # Must be > 0.
+    with pytest.raises(ValueError):
+        AlphaComplex.set_float_relative_precision(0.)
+    # Must be < 1.
+    with pytest.raises(ValueError):
+        AlphaComplex.set_float_relative_precision(1.)
+
+    points = [[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]]
+    st = AlphaComplex(points=points).create_simplex_tree()
+    filtrations = list(st.get_filtration())
+
+    # Get a better precision
+    AlphaComplex.set_float_relative_precision(1e-15)
+    assert AlphaComplex.get_float_relative_precision() == 1e-15
+
+    st = AlphaComplex(points=points).create_simplex_tree()
+    filtrations_better_resolution = list(st.get_filtration())
+
+    assert len(filtrations) == len(filtrations_better_resolution)
+    for idx in range(len(filtrations)):
+        # check simplex is the same
+        assert filtrations[idx][0] == filtrations_better_resolution[idx][0]
+        # check filtration is about the same with a relative precision of the worst case
+        assert filtrations[idx][1] == pytest.approx(filtrations_better_resolution[idx][1], rel=1e-5)
-- 
cgit v1.2.3


From f4b543bb4891b1c2aa104e45cc5a0d07b56e6977 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Thu, 31 Mar 2022 21:15:18 +0200
Subject: code review: no need to nogil for cheap functions. Repeat
 get_float_relative_precision is only for safe mode

---
 src/python/gudhi/alpha_complex.pyx | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index 9af62b43..e63796db 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -32,9 +32,9 @@ cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
         vector[double] get_point(int vertex) nogil except +
         void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except +
         @staticmethod
-        void set_float_relative_precision(double precision) nogil
+        void set_float_relative_precision(double precision)
         @staticmethod
-        double get_float_relative_precision() nogil
+        double get_float_relative_precision()
 
 # AlphaComplex python interface
 cdef class AlphaComplex:
@@ -157,7 +157,8 @@ cdef class AlphaComplex:
     def get_float_relative_precision():
         """
         :returns: The float relative precision of filtration values computation in
-            :func:`~gudhi.AlphaComplex.create_simplex_tree`.
+            :func:`~gudhi.AlphaComplex.create_simplex_tree` when the AlphaComplex is constructed with
+            :code:`precision = 'safe'` (the default).
         :rtype: float
         """
         return Alpha_complex_interface.get_float_relative_precision()
-- 
cgit v1.2.3


From 25172af7cbe5a3ea3017bed6fde53be471a7e4bf Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Tue, 5 Apr 2022 13:35:11 +0200
Subject: nogil for float_relative_precision methods

---
 src/python/gudhi/alpha_complex.pyx | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index e63796db..375e1561 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -32,9 +32,9 @@ cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
         vector[double] get_point(int vertex) nogil except +
         void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except +
         @staticmethod
-        void set_float_relative_precision(double precision)
+        void set_float_relative_precision(double precision) nogil
         @staticmethod
-        double get_float_relative_precision()
+        double get_float_relative_precision() nogil
 
 # AlphaComplex python interface
 cdef class AlphaComplex:
-- 
cgit v1.2.3


From 55b15a7c8617b770b51f61861951005a3e989327 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Wed, 6 Apr 2022 10:54:23 +0200
Subject: code review: Rewrite warnings in case of import errors

---
 src/python/gudhi/persistence_graphical_tools.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 70d550ab..3c166a56 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -97,8 +97,8 @@ def _matplotlib_can_use_tex():
         from matplotlib import checkdep_usetex
 
         return checkdep_usetex(True)
-    except ImportError:
-        warnings.warn("This function is not available, you may be missing matplotlib.")
+    except ImportError as import_error:
+        warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.")
 
 
 def plot_persistence_barcode(
@@ -233,8 +233,8 @@ def plot_persistence_barcode(
             axes.axis([axis_start, infinity, 0, ind])
         return axes
 
-    except ImportError:
-        warnings.warn("This function is not available, you may be missing matplotlib.")
+    except ImportError as import_error:
+        warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.")
 
 
 def plot_persistence_diagram(
@@ -384,8 +384,8 @@ def plot_persistence_diagram(
         axes.axis([axis_start, axis_end, axis_start, infinity + delta / 2])
         return axes
 
-    except ImportError:
-        warnings.warn("This function is not available, you may be missing matplotlib.")
+    except ImportError as import_error:
+        warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.")
 
 
 def plot_persistence_density(
@@ -549,5 +549,5 @@ def plot_persistence_density(
 
         return axes
 
-    except ImportError:
-        warnings.warn("This function is not available, you may be missing matplotlib and/or scipy.")
+    except ImportError as import_error:
+        warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.")
-- 
cgit v1.2.3


From 84a264c27e0994542f1da3a4688157568a6b69b7 Mon Sep 17 00:00:00 2001
From: Marc Glisse <marc.glisse@inria.fr>
Date: Fri, 15 Apr 2022 20:02:06 +0200
Subject: Use "and" to separate authors with bibtex

---
 biblio/how_to_cite_gudhi.bib.in | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/biblio/how_to_cite_gudhi.bib.in b/biblio/how_to_cite_gudhi.bib.in
index 54d10857..579dbf41 100644
--- a/biblio/how_to_cite_gudhi.bib.in
+++ b/biblio/how_to_cite_gudhi.bib.in
@@ -78,7 +78,7 @@
 }
 
 @incollection{gudhi:SubSampling
-, author =  "Cl\'ement Jamin, Siargey Kachanovich"
+, author =  "Cl\'ement Jamin and Siargey Kachanovich"
 , title =   "Subsampling"
 , publisher =  "{GUDHI Editorial Board}"
 , edition =     "{@GUDHI_VERSION@}"
@@ -108,7 +108,7 @@
 }
 
 @incollection{gudhi:RipsComplex
-, author =  "Cl\'ement Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse"
+, author =  "Cl\'ement Maria and Pawel Dlotko and Vincent Rouvreau and Marc Glisse"
 , title =   "Rips complex"
 , publisher =  "{GUDHI Editorial Board}"
 , edition =     "{@GUDHI_VERSION@}"
@@ -158,7 +158,7 @@
 }
 
 @incollection{gudhi:Collapse
-, author =  "Siddharth Pritam"
+, author =  "Siddharth Pritam and Marc Glisse"
 , title =   "Edge collapse"
 , publisher =  "{GUDHI Editorial Board}"
 , edition =     "{@GUDHI_VERSION@}"
-- 
cgit v1.2.3


From 5857c571388a4349934a266ca187b2c2ac10818c Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com>
Date: Thu, 21 Apr 2022 14:44:22 +0200
Subject: Make Windows compilation fails on error (#598)

Use '-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int' as a workaround
---
 .github/workflows/pip-build-windows.yml             |  8 ++++----
 .github/workflows/pip-packaging-windows.yml         |  8 ++++----
 azure-pipelines.yml                                 | 17 ++++++++++++-----
 src/cmake/modules/GUDHI_third_party_libraries.cmake |  9 +++++++++
 src/python/CMakeLists.txt                           | 12 +++++++++++-
 5 files changed, 40 insertions(+), 14 deletions(-)

diff --git a/.github/workflows/pip-build-windows.yml b/.github/workflows/pip-build-windows.yml
index 954b59d5..30b0bd94 100644
--- a/.github/workflows/pip-build-windows.yml
+++ b/.github/workflows/pip-build-windows.yml
@@ -30,14 +30,14 @@ jobs:
         run: |
           mkdir build
           cd ".\build\"
-          cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE=c:\vcpkg\scripts\buildsystems\vcpkg.cmake -DVCPKG_TARGET_TRIPLET=x64-windows ..
+          cmake -DCMAKE_BUILD_TYPE=Release -DFORCE_EIGEN_DEFAULT_DENSE_INDEX_TYPE_TO_INT=ON -DCMAKE_TOOLCHAIN_FILE=c:\vcpkg\scripts\buildsystems\vcpkg.cmake -DVCPKG_TARGET_TRIPLET=x64-windows ..
           Get-Location
           dir
           cd ".\src\python\"
-          cp "C:\vcpkg\installed\x64-windows\bin\mpfr-6.dll" ".\gudhi\"
-          cp "C:\vcpkg\installed\x64-windows\bin\gmp.dll" ".\gudhi\"
+          cp "C:\vcpkg\installed\x64-windows\bin\mpfr*.dll" ".\gudhi\"
+          cp "C:\vcpkg\installed\x64-windows\bin\gmp*.dll" ".\gudhi\"
           python setup.py bdist_wheel
-          ls dist
+          ls ".\dist\"
           cd ".\dist\"
           Get-ChildItem *.whl | ForEach-Object{python -m pip install --user $_.Name}
       - name: Test python wheel
diff --git a/.github/workflows/pip-packaging-windows.yml b/.github/workflows/pip-packaging-windows.yml
index 962ae68a..142a114c 100644
--- a/.github/workflows/pip-packaging-windows.yml
+++ b/.github/workflows/pip-packaging-windows.yml
@@ -33,14 +33,14 @@ jobs:
         run: |
           mkdir build
           cd ".\build\"
-          cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE=c:\vcpkg\scripts\buildsystems\vcpkg.cmake -DVCPKG_TARGET_TRIPLET=x64-windows ..
+          cmake -DCMAKE_BUILD_TYPE=Release -DFORCE_EIGEN_DEFAULT_DENSE_INDEX_TYPE_TO_INT=ON -DCMAKE_TOOLCHAIN_FILE=c:\vcpkg\scripts\buildsystems\vcpkg.cmake -DVCPKG_TARGET_TRIPLET=x64-windows ..
           Get-Location
           dir
           cd ".\src\python\"
-          cp "C:\vcpkg\installed\x64-windows\bin\mpfr-6.dll" ".\gudhi\"
-          cp "C:\vcpkg\installed\x64-windows\bin\gmp.dll" ".\gudhi\"
+          cp "C:\vcpkg\installed\x64-windows\bin\mpfr*.dll" ".\gudhi\"
+          cp "C:\vcpkg\installed\x64-windows\bin\gmp*.dll" ".\gudhi\"
           python setup.py bdist_wheel
-          ls dist
+          ls ".\dist\"
           cd ".\dist\"
           Get-ChildItem *.whl | ForEach-Object{python -m pip install --user $_.Name}
       - name: Test python wheel
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 21664244..a39f6d5c 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -30,7 +30,7 @@ jobs:
     - bash: |
         mkdir build
         cd build
-        cmake -DCMAKE_BUILD_TYPE:STRING=$(cmakeBuildType) -DWITH_GUDHI_TEST=ON -DWITH_GUDHI_UTILITIES=ON -DWITH_GUDHI_PYTHON=ON ..
+        cmake -DCMAKE_BUILD_TYPE:STRING=$(cmakeBuildType) -DWITH_GUDHI_EXAMPLE=ON -DWITH_GUDHI_TEST=ON -DWITH_GUDHI_UTILITIES=ON -DWITH_GUDHI_PYTHON=ON ..
         make
         make doxygen
         ctest --output-on-failure
@@ -64,20 +64,27 @@ jobs:
         vcpkg install boost-filesystem:x64-windows boost-test:x64-windows boost-program-options:x64-windows tbb:x64-windows eigen3:x64-windows cgal:x64-windows
       displayName: 'Install build dependencies'
     - script: |
-        call "C:\Program Files (x86)\Microsoft Visual Studio\2019\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" amd64
+        call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" amd64
+        IF %errorlevel% NEQ 0 exit /b %errorlevel%
         mkdir build
         cd build
-        cmake -G "Visual Studio 16 2019" -A x64 -DCMAKE_BUILD_TYPE=Release $(cmakeVcpkgFlags) $(cmakeFlags) ..
+        cmake -G "Visual Studio 17 2022" -A x64 -DCMAKE_BUILD_TYPE=Release -DFORCE_EIGEN_DEFAULT_DENSE_INDEX_TYPE_TO_INT=ON $(cmakeVcpkgFlags) $(cmakeFlags) ..
+        IF %errorlevel% NEQ 0 exit /b %errorlevel%
         MSBuild GUDHIdev.sln /m /p:Configuration=Release /p:Platform=x64
+        IF %errorlevel% NEQ 0 exit /b %errorlevel%
         ctest --output-on-failure -C Release -E diff_files
+        IF %errorlevel% NEQ 0 exit /b %errorlevel%
         cmake  -DWITH_GUDHI_PYTHON=ON .
+        IF %errorlevel% NEQ 0 exit /b %errorlevel%
         cd src\python
-        copy "C:\vcpkg\installed\x64-windows\bin\mpfr-6.dll" ".\gudhi\"
-        copy "C:\vcpkg\installed\x64-windows\bin\gmp.dll" ".\gudhi\"
+        copy "C:\vcpkg\installed\x64-windows\bin\mpfr*.dll" ".\gudhi\"
+        copy "C:\vcpkg\installed\x64-windows\bin\gmp*.dll" ".\gudhi\"
         copy "C:\vcpkg\installed\x64-windows\bin\tbb.dll" ".\gudhi\"
         copy "C:\vcpkg\installed\x64-windows\bin\tbbmalloc.dll" ".\gudhi\"
         python setup.py build_ext --inplace
+        IF %errorlevel% NEQ 0 exit /b %errorlevel%
         SET PYTHONPATH=%CD%;%PYTHONPATH%
         echo %PYTHONPATH%
         ctest --output-on-failure -C Release
+        IF %errorlevel% NEQ 0 exit /b %errorlevel%
       displayName: 'Build and test'
diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake
index 7c982b3b..6a94d1f5 100644
--- a/src/cmake/modules/GUDHI_third_party_libraries.cmake
+++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake
@@ -19,6 +19,15 @@ if(GMP_FOUND)
   endif()
 endif()
 
+# from windows vcpkg eigen 3.4.0#2 : build fails with
+# error C2440: '<function-style-cast>': cannot convert from 'Eigen::EigenBase<Derived>::Index' to '__gmp_expr<mpq_t,mpq_t>'
+# cf. https://gitlab.com/libeigen/eigen/-/issues/2476
+# Workaround is to compile with '-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int'
+if (FORCE_EIGEN_DEFAULT_DENSE_INDEX_TYPE_TO_INT)
+  message("++ User explicit demand to force EIGEN_DEFAULT_DENSE_INDEX_TYPE to int")
+  add_definitions(-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int)
+endif()
+
 # In CMakeLists.txt, when include(${CGAL_USE_FILE}), CMAKE_CXX_FLAGS are overwritten.
 # cf. http://doc.cgal.org/latest/Manual/installation.html#title40
 # A workaround is to include(${CGAL_USE_FILE}) before adding "-std=c++11".
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index fbfb0d1f..e31af02e 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -180,6 +180,15 @@ if(PYTHONINTERP_FOUND)
       set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ")
     endif ()
 
+    # from windows vcpkg eigen 3.4.0#2 : build fails with
+    # error C2440: '<function-style-cast>': cannot convert from 'Eigen::EigenBase<Derived>::Index' to '__gmp_expr<mpq_t,mpq_t>'
+    # cf. https://gitlab.com/libeigen/eigen/-/issues/2476
+    # Workaround is to compile with '-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int'
+    if (FORCE_EIGEN_DEFAULT_DENSE_INDEX_TYPE_TO_INT)
+      set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int', ")
+    endif()
+
+
     add_gudhi_debug_info("Boost version ${Boost_VERSION}")
     if(CGAL_FOUND)
       if(NOT CGAL_VERSION VERSION_LESS 5.3.0)
@@ -215,13 +224,14 @@ if(PYTHONINTERP_FOUND)
         endif(NOT GMP_LIBRARIES_DIR)
         add_GUDHI_PYTHON_lib_dir(${GMP_LIBRARIES_DIR})
         message("** Add gmp ${GMP_LIBRARIES_DIR}")
+        # When FORCE_CGAL_NOT_TO_BUILD_WITH_GMPXX is set, not defining CGAL_USE_GMPXX is sufficient enough
         if(GMPXX_FOUND)
           add_gudhi_debug_info("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}")
           set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMPXX', ")
           add_GUDHI_PYTHON_lib("${GMPXX_LIBRARIES}")
           add_GUDHI_PYTHON_lib_dir(${GMPXX_LIBRARIES_DIR})
           message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}")
-        endif(GMPXX_FOUND)
+        endif()
       endif(GMP_FOUND)
       if(MPFR_FOUND)
         add_gudhi_debug_info("MPFR_LIBRARIES = ${MPFR_LIBRARIES}")
-- 
cgit v1.2.3


From 5e9811c4d79436f18868a2f2128ef5d0310e66ce Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Wed, 27 Apr 2022 09:41:52 +0200
Subject: latest

---
 azure-pipelines.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index a39f6d5c..31264c37 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -5,7 +5,7 @@ jobs:
     timeoutInMinutes: 0
     cancelTimeoutInMinutes: 60
     pool:
-      vmImage: macOS-10.15
+      vmImage: macOS-latest
     variables:
       pythonVersion: '3.7'
       cmakeBuildType: Release
-- 
cgit v1.2.3


From 79278cfc0aecbfdbf8a66e79ccef44534abb2399 Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Fri, 29 Apr 2022 17:48:19 +0200
Subject: plot_persistence_diagram and plot_persistence_barcode improvements

---
 src/python/gudhi/persistence_graphical_tools.py | 68 ++++++-------------------
 1 file changed, 16 insertions(+), 52 deletions(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 3c166a56..604018d1 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -105,7 +105,7 @@ def plot_persistence_barcode(
     persistence=[],
     persistence_file="",
     alpha=0.6,
-    max_intervals=1000,
+    max_intervals=20000,
     inf_delta=0.1,
     legend=False,
     colormap=None,
@@ -127,7 +127,7 @@ def plot_persistence_barcode(
     :type alpha: float.
     :param max_intervals: maximal number of intervals to display.
         Selected intervals are those with the longest life time. Set it
-        to 0 to see all. Default value is 1000.
+        to 0 to see all. Default value is 20000.
     :type max_intervals: int.
     :param inf_delta: Infinity is placed at :code:`((max_death - min_birth) x
         inf_delta)` above :code:`max_death` value. A reasonable value is
@@ -189,36 +189,11 @@ def plot_persistence_barcode(
             _, axes = plt.subplots(1, 1)
         if colormap == None:
             colormap = plt.cm.Set1.colors
-        ind = 0
-
-        # Draw horizontal bars in loop
-        for interval in reversed(persistence):
-            try:
-                if float(interval[1][1]) != float("inf"):
-                    # Finite death case
-                    axes.barh(
-                        ind,
-                        (interval[1][1] - interval[1][0]),
-                        height=0.8,
-                        left=interval[1][0],
-                        alpha=alpha,
-                        color=colormap[interval[0]],
-                        linewidth=0,
-                    )
-                else:
-                    # Infinite death case for diagram to be nicer
-                    axes.barh(
-                        ind,
-                        (infinity - interval[1][0]),
-                        height=0.8,
-                        left=interval[1][0],
-                        alpha=alpha,
-                        color=colormap[interval[0]],
-                        linewidth=0,
-                    )
-                ind = ind + 1
-            except IndexError:
-                pass
+
+        x=[birth for (dim,(birth,death)) in persistence]
+        y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence]
+        c=[colormap[dim] for (dim,(birth,death)) in persistence]
+        axes.barh(list(reversed(range(len(x)))), y, height=0.8, left=x, alpha=alpha, color=c, linewidth=0)
 
         if legend:
             dimensions = list(set(item[0] for item in persistence))
@@ -229,8 +204,8 @@ def plot_persistence_barcode(
         axes.set_title("Persistence barcode", fontsize=fontsize)
 
         # Ends plot on infinity value and starts a little bit before min_birth
-        if ind != 0:
-            axes.axis([axis_start, infinity, 0, ind])
+        if len(x) != 0:
+            axes.axis([axis_start, infinity, 0, len(x)])
         return axes
 
     except ImportError as import_error:
@@ -242,7 +217,7 @@ def plot_persistence_diagram(
     persistence_file="",
     alpha=0.6,
     band=0.0,
-    max_intervals=1000,
+    max_intervals=1000000,
     inf_delta=0.1,
     legend=False,
     colormap=None,
@@ -266,7 +241,7 @@ def plot_persistence_diagram(
     :type band: float.
     :param max_intervals: maximal number of intervals to display.
         Selected intervals are those with the longest life time. Set it
-        to 0 to see all. Default value is 1000.
+        to 0 to see all. Default value is 1000000.
     :type max_intervals: int.
     :param inf_delta: Infinity is placed at :code:`((max_death - min_birth) x
         inf_delta)` above :code:`max_death` value. A reasonable value is
@@ -346,22 +321,11 @@ def plot_persistence_diagram(
         # line display of equation : birth = death
         axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
 
-        # Draw points in loop
-        pts_at_infty = False  # Records presence of pts at infty
-        for interval in reversed(persistence):
-            try:
-                if float(interval[1][1]) != float("inf"):
-                    # Finite death case
-                    axes.scatter(
-                        interval[1][0], interval[1][1], alpha=alpha, color=colormap[interval[0]],
-                    )
-                else:
-                    pts_at_infty = True
-                    # Infinite death case for diagram to be nicer
-                    axes.scatter(interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]])
-            except IndexError:
-                pass
-        if pts_at_infty:
+        x=[birth for (dim,(birth,death)) in persistence]
+        y=[death if death != float("inf") else infinity for (dim,(birth,death)) in persistence]
+        c=[colormap[dim] for (dim,(birth,death)) in persistence]
+        axes.scatter(x,y,alpha=alpha,color=c)
+        if float("inf") in (death for (dim,(birth,death)) in persistence):
             # infinity line and text
             axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha)
             # Infinity label
-- 
cgit v1.2.3


From 3952b61604db976f4864b85ebdafa0962766b8bc Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Tue, 10 May 2022 10:00:02 +0200
Subject: Fix limit tests for plot (and warning in test)

---
 src/python/gudhi/persistence_graphical_tools.py    | 32 +++++++++++++++-------
 .../test/test_persistence_graphical_tools.py       |  6 +++-
 2 files changed, 27 insertions(+), 11 deletions(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 604018d1..930df825 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -190,12 +190,17 @@ def plot_persistence_barcode(
         if colormap == None:
             colormap = plt.cm.Set1.colors
 
-        x=[birth for (dim,(birth,death)) in persistence]
-        y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence]
-        c=[colormap[dim] for (dim,(birth,death)) in persistence]
+        non_empty_diagram = len(persistence[0]) > 0
+        if non_empty_diagram:
+            x=[birth for (dim,(birth,death)) in persistence]
+            y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence]
+            c=[colormap[dim] for (dim,(birth,death)) in persistence]
+        else:
+            x, y, c = [], [], []
+
         axes.barh(list(reversed(range(len(x)))), y, height=0.8, left=x, alpha=alpha, color=c, linewidth=0)
 
-        if legend:
+        if non_empty_diagram and legend:
             dimensions = list(set(item[0] for item in persistence))
             axes.legend(
                 handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions], loc="lower right",
@@ -321,11 +326,16 @@ def plot_persistence_diagram(
         # line display of equation : birth = death
         axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
 
-        x=[birth for (dim,(birth,death)) in persistence]
-        y=[death if death != float("inf") else infinity for (dim,(birth,death)) in persistence]
-        c=[colormap[dim] for (dim,(birth,death)) in persistence]
+        non_empty_diagram = len(persistence[0]) > 0
+        if non_empty_diagram:
+            x=[birth for (dim,(birth,death)) in persistence]
+            y=[death if death != float("inf") else infinity for (dim,(birth,death)) in persistence]
+            c=[colormap[dim] for (dim,(birth,death)) in persistence]
+        else:
+            x, y, c = [], [], []
+
         axes.scatter(x,y,alpha=alpha,color=c)
-        if float("inf") in (death for (dim,(birth,death)) in persistence):
+        if non_empty_diagram and float("inf") in (death for (dim,(birth,death)) in persistence):
             # infinity line and text
             axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha)
             # Infinity label
@@ -337,7 +347,7 @@ def plot_persistence_diagram(
             axes.set_yticks(yt)
             axes.set_yticklabels(ytl)
 
-        if legend:
+        if non_empty_diagram and legend:
             dimensions = list(set(item[0] for item in persistence))
             axes.legend(handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions])
 
@@ -482,6 +492,7 @@ def plot_persistence_density(
             zi = k(np.vstack([xi.flatten(), yi.flatten()]))
             # Make the plot
             img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading="auto")
+            non_empty_diagram = True
 
         # IndexError on empty diagrams, ValueError on only inf death values
         except (IndexError, ValueError):
@@ -489,6 +500,7 @@ def plot_persistence_density(
             birth_max = 1.0
             death_min = 0.0
             death_max = 1.0
+            non_empty_diagram = False
             pass
 
         # line display of equation : birth = death
@@ -504,7 +516,7 @@ def plot_persistence_density(
                 )
             )
 
-        if legend:
+        if non_empty_diagram and legend:
             plt.colorbar(img, ax=axes)
 
         axes.set_xlabel("Birth", fontsize=fontsize)
diff --git a/src/python/test/test_persistence_graphical_tools.py b/src/python/test/test_persistence_graphical_tools.py
index 7d9bae90..994791d7 100644
--- a/src/python/test/test_persistence_graphical_tools.py
+++ b/src/python/test/test_persistence_graphical_tools.py
@@ -15,7 +15,7 @@ import pytest
 
 
 def test_array_handler():
-    diags = np.array([[1, 2], [3, 4], [5, 6]], np.float)
+    diags = np.array([[1, 2], [3, 4], [5, 6]], float)
     arr_diags = gd.persistence_graphical_tools._array_handler(diags)
     for idx in range(len(diags)):
         assert arr_diags[idx][0] == 0
@@ -98,8 +98,12 @@ def test_limit_to_max_intervals():
 def _limit_plot_persistence(function):
     pplot = function(persistence=[()])
     assert issubclass(type(pplot), plt.axes.SubplotBase)
+    pplot = function(persistence=[()], legend=True)
+    assert issubclass(type(pplot), plt.axes.SubplotBase)
     pplot = function(persistence=[(0, float("inf"))])
     assert issubclass(type(pplot), plt.axes.SubplotBase)
+    pplot = function(persistence=[(0, float("inf"))], legend=True)
+    assert issubclass(type(pplot), plt.axes.SubplotBase)
 
 
 def test_limit_plot_persistence():
-- 
cgit v1.2.3


From 393216cc99af79526b1bb90ae4da21116884bcbd Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 16 May 2022 13:47:16 +0200
Subject: code review: limit test was not respecting the one described in the
 documentation

---
 src/python/gudhi/persistence_graphical_tools.py    | 51 +++++++++-------------
 .../test/test_persistence_graphical_tools.py       |  4 +-
 2 files changed, 23 insertions(+), 32 deletions(-)

diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 930df825..7ed11360 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -190,17 +190,13 @@ def plot_persistence_barcode(
         if colormap == None:
             colormap = plt.cm.Set1.colors
 
-        non_empty_diagram = len(persistence[0]) > 0
-        if non_empty_diagram:
-            x=[birth for (dim,(birth,death)) in persistence]
-            y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence]
-            c=[colormap[dim] for (dim,(birth,death)) in persistence]
-        else:
-            x, y, c = [], [], []
+        x=[birth for (dim,(birth,death)) in persistence]
+        y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence]
+        c=[colormap[dim] for (dim,(birth,death)) in persistence]
 
         axes.barh(list(reversed(range(len(x)))), y, height=0.8, left=x, alpha=alpha, color=c, linewidth=0)
 
-        if non_empty_diagram and legend:
+        if legend:
             dimensions = list(set(item[0] for item in persistence))
             axes.legend(
                 handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions], loc="lower right",
@@ -326,16 +322,12 @@ def plot_persistence_diagram(
         # line display of equation : birth = death
         axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
 
-        non_empty_diagram = len(persistence[0]) > 0
-        if non_empty_diagram:
-            x=[birth for (dim,(birth,death)) in persistence]
-            y=[death if death != float("inf") else infinity for (dim,(birth,death)) in persistence]
-            c=[colormap[dim] for (dim,(birth,death)) in persistence]
-        else:
-            x, y, c = [], [], []
+        x=[birth for (dim,(birth,death)) in persistence]
+        y=[death if death != float("inf") else infinity for (dim,(birth,death)) in persistence]
+        c=[colormap[dim] for (dim,(birth,death)) in persistence]
 
         axes.scatter(x,y,alpha=alpha,color=c)
-        if non_empty_diagram and float("inf") in (death for (dim,(birth,death)) in persistence):
+        if float("inf") in (death for (dim,(birth,death)) in persistence):
             # infinity line and text
             axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha)
             # Infinity label
@@ -347,7 +339,7 @@ def plot_persistence_diagram(
             axes.set_yticks(yt)
             axes.set_yticklabels(ytl)
 
-        if non_empty_diagram and legend:
+        if legend:
             dimensions = list(set(item[0] for item in persistence))
             axes.legend(handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions])
 
@@ -459,16 +451,15 @@ def plot_persistence_density(
             _, axes = plt.subplots(1, 1)
 
         try:
-            if len(persistence) > 0:
-                # if not read from file but given by an argument
-                persistence = _array_handler(persistence)
-                persistence_dim = np.array(
-                    [
-                        (dim_interval[1][0], dim_interval[1][1])
-                        for dim_interval in persistence
-                        if (dim_interval[0] == dimension) or (dimension is None)
-                    ]
-                )
+            # if not read from file but given by an argument
+            persistence = _array_handler(persistence)
+            persistence_dim = np.array(
+                [
+                    (dim_interval[1][0], dim_interval[1][1])
+                    for dim_interval in persistence
+                    if (dim_interval[0] == dimension) or (dimension is None)
+                ]
+            )
             persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])]
             persistence_dim = np.array(
                 _limit_to_max_intervals(
@@ -492,7 +483,7 @@ def plot_persistence_density(
             zi = k(np.vstack([xi.flatten(), yi.flatten()]))
             # Make the plot
             img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading="auto")
-            non_empty_diagram = True
+            plot_success = True
 
         # IndexError on empty diagrams, ValueError on only inf death values
         except (IndexError, ValueError):
@@ -500,7 +491,7 @@ def plot_persistence_density(
             birth_max = 1.0
             death_min = 0.0
             death_max = 1.0
-            non_empty_diagram = False
+            plot_success = False
             pass
 
         # line display of equation : birth = death
@@ -516,7 +507,7 @@ def plot_persistence_density(
                 )
             )
 
-        if non_empty_diagram and legend:
+        if plot_success and legend:
             plt.colorbar(img, ax=axes)
 
         axes.set_xlabel("Birth", fontsize=fontsize)
diff --git a/src/python/test/test_persistence_graphical_tools.py b/src/python/test/test_persistence_graphical_tools.py
index 994791d7..44decdd7 100644
--- a/src/python/test/test_persistence_graphical_tools.py
+++ b/src/python/test/test_persistence_graphical_tools.py
@@ -96,9 +96,9 @@ def test_limit_to_max_intervals():
 
 
 def _limit_plot_persistence(function):
-    pplot = function(persistence=[()])
+    pplot = function(persistence=[])
     assert issubclass(type(pplot), plt.axes.SubplotBase)
-    pplot = function(persistence=[()], legend=True)
+    pplot = function(persistence=[], legend=True)
     assert issubclass(type(pplot), plt.axes.SubplotBase)
     pplot = function(persistence=[(0, float("inf"))])
     assert issubclass(type(pplot), plt.axes.SubplotBase)
-- 
cgit v1.2.3


From 563295694afda3dfa37cccb63a73a9f22131480e Mon Sep 17 00:00:00 2001
From: Vincent Rouvreau <vincent.rouvreau@inria.fr>
Date: Mon, 16 May 2022 17:14:29 +0200
Subject: code review: use isinstance(x, T) instead of issubclass(type(x),T)

---
 src/python/test/test_persistence_graphical_tools.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/python/test/test_persistence_graphical_tools.py b/src/python/test/test_persistence_graphical_tools.py
index 44decdd7..c19836b7 100644
--- a/src/python/test/test_persistence_graphical_tools.py
+++ b/src/python/test/test_persistence_graphical_tools.py
@@ -97,13 +97,13 @@ def test_limit_to_max_intervals():
 
 def _limit_plot_persistence(function):
     pplot = function(persistence=[])
-    assert issubclass(type(pplot), plt.axes.SubplotBase)
+    assert isinstance(pplot, plt.axes.SubplotBase)
     pplot = function(persistence=[], legend=True)
-    assert issubclass(type(pplot), plt.axes.SubplotBase)
+    assert isinstance(pplot, plt.axes.SubplotBase)
     pplot = function(persistence=[(0, float("inf"))])
-    assert issubclass(type(pplot), plt.axes.SubplotBase)
+    assert isinstance(pplot, plt.axes.SubplotBase)
     pplot = function(persistence=[(0, float("inf"))], legend=True)
-    assert issubclass(type(pplot), plt.axes.SubplotBase)
+    assert isinstance(pplot, plt.axes.SubplotBase)
 
 
 def test_limit_plot_persistence():
-- 
cgit v1.2.3