Merge pull request #331 from mglisse/insert

convert array of filtration values to SimplexTree

6 files changed, 359 insertions, 151 deletions

diff --git a/.github/next_release.md b/.github/next_release.md
index d5fcef1c..929a7ce6 100644
--- a/.github/next_release.md
+++ b/.github/next_release.md
@@ -9,6 +9,9 @@ Below is a list of changes made since GUDHI 3.6.0:
 - [Module](link)
      - ...
 
+- [Simplex tree](https://gudhi.inria.fr/python/latest/simplex_tree_ref.html)
+     - New functions to initialize from a matrix or insert batches of simplices of the same dimension.
+
 - [Rips complex](https://gudhi.inria.fr/python/latest/rips_complex_user.html)
      - Construction now rejects positional arguments, you need to specify `points=X`.
 
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 9059219c..ef9f8428 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -24,6 +24,7 @@
 #include <boost/iterator/transform_iterator.hpp>
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/range/adaptor/reversed.hpp>
+#include <boost/range/size.hpp>
 #include <boost/container/static_vector.hpp>
 
 #ifdef GUDHI_USE_TBB
@@ -702,10 +703,10 @@ class Simplex_tree {
     return true;
   }
 
- private:
-  /** \brief Inserts a simplex represented by a vector of vertex.
-   * @param[in]  simplex    vector of Vertex_handles, representing the vertices of the new simplex. The vector must be
-   * sorted by increasing vertex handle order.
+ protected:
+  /** \brief Inserts a simplex represented by a range of vertex.
+   * @param[in]  simplex    range of Vertex_handles, representing the vertices of the new simplex. The range must be
+   * sorted by increasing vertex handle order, and not empty.
    * @param[in]  filtration the filtration value assigned to the new simplex.
    * @return If the new simplex is inserted successfully (i.e. it was not in the
    * simplicial complex yet) the bool is set to true and the Simplex_handle is the handle assigned
@@ -717,12 +718,13 @@ class Simplex_tree {
    * null_simplex.
    * 
   */
-  std::pair<Simplex_handle, bool> insert_vertex_vector(const std::vector<Vertex_handle>& simplex,
+  template <class RandomVertexHandleRange = std::initializer_list<Vertex_handle>>
+  std::pair<Simplex_handle, bool> insert_simplex_raw(const RandomVertexHandleRange& simplex,
                                                      Filtration_value filtration) {
     Siblings * curr_sib = &root_;
     std::pair<Simplex_handle, bool> res_insert;
     auto vi = simplex.begin();
-    for (; vi != simplex.end() - 1; ++vi) {
+    for (; vi != std::prev(simplex.end()); ++vi) {
       GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex");
       res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
       if (!(has_children(res_insert.first))) {
@@ -743,9 +745,10 @@ class Simplex_tree {
       return std::pair<Simplex_handle, bool>(null_simplex(), false);
     }
     // otherwise the insertion has succeeded - size is a size_type
-    if (static_cast<int>(simplex.size()) - 1 > dimension_) {
+    int dim = static_cast<int>(boost::size(simplex)) - 1;
+    if (dim > dimension_) {
       // Update dimension if needed
-      dimension_ = static_cast<int>(simplex.size()) - 1;
+      dimension_ = dim;
     }
     return res_insert;
   }
@@ -786,7 +789,7 @@ class Simplex_tree {
     // Copy before sorting
     std::vector<Vertex_handle> copy(first, last);
     std::sort(std::begin(copy), std::end(copy));
-    return insert_vertex_vector(copy, filtration);
+    return insert_simplex_raw(copy, filtration);
   }
 
     /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of
@@ -1598,7 +1601,7 @@ class Simplex_tree {
     Simplex_tree st_copy = *this;
 
     // Add point for coning the simplicial complex
-    this->insert_simplex({maxvert}, -3);
+    this->insert_simplex_raw({maxvert}, -3);
 
     // For each simplex
     std::vector<Vertex_handle> vr;
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 5642f82d..f86f1232 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -56,6 +56,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         int upper_bound_dimension() nogil
         bool find_simplex(vector[int] simplex) nogil
         bool insert(vector[int] simplex, double filtration) nogil
+        void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) nogil
         vector[pair[vector[int], double]] get_star(vector[int] simplex) nogil
         vector[pair[vector[int], double]] get_cofaces(vector[int] simplex, int dimension) nogil
         void expansion(int max_dim) nogil except +
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 05bfe22e..24b970c4 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -8,14 +8,23 @@
 #   - YYYY/MM Author: Description of the modification
 
 from cython.operator import dereference, preincrement
-from libc.stdint cimport intptr_t
+from libc.stdint cimport intptr_t, int32_t, int64_t
 import numpy as np
 cimport gudhi.simplex_tree
+cimport cython
 
 __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
 __license__ = "MIT"
 
+ctypedef fused some_int:
+    int32_t
+    int64_t
+
+ctypedef fused some_float:
+    float
+    double
+
 cdef bool callback(vector[int] simplex, void *blocker_func):
     return (<object>blocker_func)(simplex)
 
@@ -228,6 +237,87 @@ cdef class SimplexTree:
         """
         return self.get_ptr().insert(simplex, <double>filtration)
 
+    @staticmethod
+    @cython.boundscheck(False)
+    def create_from_array(filtrations, double max_filtration=np.inf):
+        """Creates a new, empty complex and inserts vertices and edges. The vertices are numbered from 0 to n-1, and
+        the filtration values are encoded in the array, with the diagonal representing the vertices. It is the
+        caller's responsibility to ensure that this defines a filtration, which can be achieved with either::
+
+            filtrations[np.diag_indices_from(filtrations)] = filtrations.min(axis=1)
+
+        or::
+
+            diag = filtrations.diagonal()
+            filtrations = np.fmax(np.fmax(filtrations, diag[:, None]), diag[None, :])
+
+        :param filtrations: the filtration values of the vertices and edges to insert. The matrix is assumed to be symmetric.
+        :type filtrations: numpy.ndarray of shape (n,n)
+        :param max_filtration: only insert vertices and edges with filtration values no larger than max_filtration
+        :type max_filtration: float
+        :returns: the new complex
+        :rtype: SimplexTree
+        """
+        # TODO: document which half of the matrix is actually read?
+        filtrations = np.asanyarray(filtrations, dtype=float)
+        cdef double[:,:] F = filtrations
+        ret = SimplexTree()
+        cdef int n = F.shape[0]
+        assert n == F.shape[1], 'create_from_array() expects a square array'
+        with nogil:
+            ret.get_ptr().insert_matrix(&F[0,0], n, F.strides[0], F.strides[1], max_filtration)
+        return ret
+
+    def insert_edges_from_coo_matrix(self, edges):
+        """Inserts edges given by a sparse matrix in `COOrdinate format
+        <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html>`_.
+        If an edge is repeated, the smallest filtration value is used. Missing entries are not inserted.
+        Diagonal entries are currently interpreted as vertices, although we do not guarantee this behavior
+        in the future, and this is only useful if you want to insert vertices with a smaller filtration value
+        than the smallest edge containing it, since vertices are implicitly inserted together with the edges.
+
+        :param edges: the edges to insert and their filtration values.
+        :type edges: scipy.sparse.coo_matrix of shape (n,n)
+
+        .. seealso:: :func:`insert_batch`
+        """
+        # TODO: optimize this?
+        for edge in zip(edges.row, edges.col, edges.data):
+            self.get_ptr().insert((edge[0], edge[1]), edge[2])
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def insert_batch(self, some_int[:,:] vertex_array, some_float[:] filtrations):
+        """Inserts k-simplices given by a sparse array in a format similar
+        to `torch.sparse <https://pytorch.org/docs/stable/sparse.html>`_.
+        The n-th simplex has vertices `vertex_array[0,n]`, ...,
+        `vertex_array[k,n]` and filtration value `filtrations[n]`.
+        If a simplex is repeated, the smallest filtration value is used.
+        Simplices with a repeated vertex are currently interpreted as lower
+        dimensional simplices, but we do not guarantee this behavior in the
+        future. Any time a simplex is inserted, its faces are inserted as well
+        if needed to preserve a simplicial complex.
+
+        :param vertex_array: the k-simplices to insert.
+        :type vertex_array: numpy.array of shape (k+1,n)
+        :param filtrations: the filtration values.
+        :type filtrations: numpy.array of shape (n,)
+        """
+        # This may be slow if we end up inserting vertices in a bad order (flat_map).
+        # We could first insert the vertices from np.unique(vertex_array), or leave it to the caller.
+        cdef Py_ssize_t k = vertex_array.shape[0]
+        cdef Py_ssize_t n = vertex_array.shape[1]
+        assert filtrations.shape[0] == n, 'inconsistent sizes for vertex_array and filtrations'
+        cdef Py_ssize_t i
+        cdef Py_ssize_t j
+        cdef vector[int] v
+        with nogil:
+            for i in range(n):
+                for j in range(k):
+                    v.push_back(vertex_array[j, i])
+                self.get_ptr().insert(v, filtrations[i])
+                v.clear()
+
     def get_simplices(self):
         """This function returns a generator with simplices and their given
         filtration values.
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index 3848c5ad..0317ea39 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -40,6 +40,8 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
   using Complex_simplex_iterator = typename Base::Complex_simplex_iterator;
   using Extended_filtration_data = typename Base::Extended_filtration_data;
   using Boundary_simplex_iterator = typename Base::Boundary_simplex_iterator;
+  using Siblings = typename Base::Siblings;
+  using Node = typename Base::Node;
   typedef bool (*blocker_func_t)(Simplex simplex, void *user_data);
 
  public:
@@ -62,6 +64,30 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     return (result.second);
   }
 
+  void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) {
+    // We could delegate to insert_graph, but wrapping the matrix in a graph interface is too much work,
+    // and this is a bit more efficient.
+    auto& rm = this->root()->members_;
+    for(int i=0; i<n; ++i) {
+      char* p = reinterpret_cast<char*>(filtrations) + i * stride0;
+      double fv = *reinterpret_cast<double*>(p + i * stride1);
+      if(fv > max_filtration) continue;
+      auto sh = rm.emplace_hint(rm.end(), i, Node(this->root(), fv));
+      Siblings* children = nullptr;
+      // Should we make a first pass to count the number of edges so we can reserve the right space?
+      for(int j=i+1; j<n; ++j) {
+        double fe = *reinterpret_cast<double*>(p + j * stride1);
+        if(fe > max_filtration) continue;
+        if(!children) {
+          children = new Siblings(this->root(), i);
+          sh->second.assign_children(children);
+        }
+        children->members().emplace_hint(children->members().end(), j, Node(children, fe));
+      }
+    }
+
+  }
+
   // Do not interface this function, only used in alpha complex interface for complex creation
   bool insert_simplex(const Simplex& simplex, Filtration_value filtration = 0) {
     Insertion_result result = Base::insert_simplex(simplex, filtration);
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 54bafed5..2ccbfbf5 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -249,6 +249,7 @@ def test_make_filtration_non_decreasing():
     assert st.filtration([3, 4]) == 2.0
     assert st.filtration([4, 5]) == 2.0
 
+
 def test_extend_filtration():
 
     # Inserted simplex:
@@ -257,86 +258,87 @@ def test_extend_filtration():
     #     / \ /
     #    o   o
     #   /2\ /3
-    #  o   o        
-    #  1   0        
-
-    st = SimplexTree()                                                                                                                     
-    st.insert([0,2])
-    st.insert([1,2])
-    st.insert([0,3])
-    st.insert([2,5])
-    st.insert([3,4])
-    st.insert([3,5])                                                                                                         
-    st.assign_filtration([0], 1.)                                                                                                                
-    st.assign_filtration([1], 2.)                                                                                                                
-    st.assign_filtration([2], 3.)                                                                                                                
-    st.assign_filtration([3], 4.)                                                                                                                
-    st.assign_filtration([4], 5.)                                                                                                                
-    st.assign_filtration([5], 6.)                                                                                                                
-
-    assert list(st.get_filtration()) == [                                                                                                                         
-        ([0, 2], 0.0), 
-        ([1, 2], 0.0), 
-        ([0, 3], 0.0), 
-        ([3, 4], 0.0), 
-        ([2, 5], 0.0), 
-        ([3, 5], 0.0), 
-        ([0], 1.0), 
-        ([1], 2.0), 
-        ([2], 3.0), 
-        ([3], 4.0), 
-        ([4], 5.0), 
-        ([5], 6.0)
+    #  o   o
+    #  1   0
+
+    st = SimplexTree()
+    st.insert([0, 2])
+    st.insert([1, 2])
+    st.insert([0, 3])
+    st.insert([2, 5])
+    st.insert([3, 4])
+    st.insert([3, 5])
+    st.assign_filtration([0], 1.0)
+    st.assign_filtration([1], 2.0)
+    st.assign_filtration([2], 3.0)
+    st.assign_filtration([3], 4.0)
+    st.assign_filtration([4], 5.0)
+    st.assign_filtration([5], 6.0)
+
+    assert list(st.get_filtration()) == [
+        ([0, 2], 0.0),
+        ([1, 2], 0.0),
+        ([0, 3], 0.0),
+        ([3, 4], 0.0),
+        ([2, 5], 0.0),
+        ([3, 5], 0.0),
+        ([0], 1.0),
+        ([1], 2.0),
+        ([2], 3.0),
+        ([3], 4.0),
+        ([4], 5.0),
+        ([5], 6.0),
     ]
-        
+
     st.extend_filtration()
-    
-    assert list(st.get_filtration()) == [                                                                                                                         
-        ([6], -3.0), 
-        ([0], -2.0), 
-        ([1], -1.8), 
-        ([2], -1.6), 
-        ([0, 2], -1.6), 
-        ([1, 2], -1.6), 
-        ([3], -1.4), 
-        ([0, 3], -1.4), 
-        ([4], -1.2), 
-        ([3, 4], -1.2), 
-        ([5], -1.0), 
-        ([2, 5], -1.0), 
-        ([3, 5], -1.0), 
-        ([5, 6], 1.0), 
-        ([4, 6], 1.2), 
-        ([3, 6], 1.4), 
+
+    assert list(st.get_filtration()) == [
+        ([6], -3.0),
+        ([0], -2.0),
+        ([1], -1.8),
+        ([2], -1.6),
+        ([0, 2], -1.6),
+        ([1, 2], -1.6),
+        ([3], -1.4),
+        ([0, 3], -1.4),
+        ([4], -1.2),
+        ([3, 4], -1.2),
+        ([5], -1.0),
+        ([2, 5], -1.0),
+        ([3, 5], -1.0),
+        ([5, 6], 1.0),
+        ([4, 6], 1.2),
+        ([3, 6], 1.4),
         ([3, 4, 6], 1.4),
-        ([3, 5, 6], 1.4), 
-        ([2, 6], 1.6), 
-        ([2, 5, 6], 1.6), 
-        ([1, 6], 1.8), 
-        ([1, 2, 6], 1.8), 
-        ([0, 6], 2.0), 
-        ([0, 2, 6], 2.0), 
-        ([0, 3, 6], 2.0)
+        ([3, 5, 6], 1.4),
+        ([2, 6], 1.6),
+        ([2, 5, 6], 1.6),
+        ([1, 6], 1.8),
+        ([1, 2, 6], 1.8),
+        ([0, 6], 2.0),
+        ([0, 2, 6], 2.0),
+        ([0, 3, 6], 2.0),
     ]
 
-    dgms = st.extended_persistence(min_persistence=-1.)
+    dgms = st.extended_persistence(min_persistence=-1.0)
     assert len(dgms) == 4
     # Sort by (death-birth) descending - we are only interested in those with the longest life span
     for idx in range(4):
-        dgms[idx] = sorted(dgms[idx], key=lambda x:(-abs(x[1][0]-x[1][1])))
+        dgms[idx] = sorted(dgms[idx], key=lambda x: (-abs(x[1][0] - x[1][1])))
+
+    assert dgms[0][0][1][0] == pytest.approx(2.0)
+    assert dgms[0][0][1][1] == pytest.approx(3.0)
+    assert dgms[1][0][1][0] == pytest.approx(5.0)
+    assert dgms[1][0][1][1] == pytest.approx(4.0)
+    assert dgms[2][0][1][0] == pytest.approx(1.0)
+    assert dgms[2][0][1][1] == pytest.approx(6.0)
+    assert dgms[3][0][1][0] == pytest.approx(6.0)
+    assert dgms[3][0][1][1] == pytest.approx(1.0)
 
-    assert dgms[0][0][1][0] == pytest.approx(2.)
-    assert dgms[0][0][1][1] == pytest.approx(3.)
-    assert dgms[1][0][1][0] == pytest.approx(5.)
-    assert dgms[1][0][1][1] == pytest.approx(4.)
-    assert dgms[2][0][1][0] == pytest.approx(1.)
-    assert dgms[2][0][1][1] == pytest.approx(6.)
-    assert dgms[3][0][1][0] == pytest.approx(6.)
-    assert dgms[3][0][1][1] == pytest.approx(1.) 
 
 def test_simplices_iterator():
     st = SimplexTree()
-    
+
     assert st.insert([0, 1, 2], filtration=4.0) == True
     assert st.insert([2, 3, 4], filtration=2.0) == True
 
@@ -346,9 +348,10 @@ def test_simplices_iterator():
         print("filtration is: ", simplex[1])
         assert st.filtration(simplex[0]) == simplex[1]
 
+
 def test_collapse_edges():
     st = SimplexTree()
-    
+
     assert st.insert([0, 1], filtration=1.0) == True
     assert st.insert([1, 2], filtration=1.0) == True
     assert st.insert([2, 3], filtration=1.0) == True
@@ -360,31 +363,33 @@ def test_collapse_edges():
 
     st.collapse_edges()
     assert st.num_simplices() == 9
-    assert st.find([0, 2]) == False # [1, 3] would be fine as well
+    assert st.find([0, 2]) == False  # [1, 3] would be fine as well
     for simplex in st.get_skeleton(0):
-        assert simplex[1] == 1.
+        assert simplex[1] == 1.0
+
 
 def test_reset_filtration():
     st = SimplexTree()
-    
-    assert st.insert([0, 1, 2], 3.) == True
-    assert st.insert([0, 3], 2.) == True
-    assert st.insert([3, 4, 5], 3.) == True
-    assert st.insert([0, 1, 6, 7], 4.) == True
+
+    assert st.insert([0, 1, 2], 3.0) == True
+    assert st.insert([0, 3], 2.0) == True
+    assert st.insert([3, 4, 5], 3.0) == True
+    assert st.insert([0, 1, 6, 7], 4.0) == True
 
     # Guaranteed by construction
     for simplex in st.get_simplices():
-        assert st.filtration(simplex[0]) >= 2.
-    
+        assert st.filtration(simplex[0]) >= 2.0
+
     # dimension until 5 even if simplex tree is of dimension 3 to test the limits
     for dimension in range(5, -1, -1):
-        st.reset_filtration(0., dimension)
+        st.reset_filtration(0.0, dimension)
         for simplex in st.get_skeleton(3):
             print(simplex)
             if len(simplex[0]) < (dimension) + 1:
-                assert st.filtration(simplex[0]) >= 2.
+                assert st.filtration(simplex[0]) >= 2.0
             else:
-                assert st.filtration(simplex[0]) == 0.
+                assert st.filtration(simplex[0]) == 0.0
+
 
 def test_boundaries_iterator():
     st = SimplexTree()
@@ -400,16 +405,17 @@ def test_boundaries_iterator():
         list(st.get_boundaries([]))
 
     with pytest.raises(RuntimeError):
-        list(st.get_boundaries([0, 4])) # (0, 4) does not exist
+        list(st.get_boundaries([0, 4]))  # (0, 4) does not exist
 
     with pytest.raises(RuntimeError):
-        list(st.get_boundaries([6])) # (6) does not exist
+        list(st.get_boundaries([6]))  # (6) does not exist
+
 
 def test_persistence_intervals_in_dimension():
     # Here is our triangulation of a 2-torus - taken from https://dioscuri-tda.org/Paris_TDA_Tutorial_2021.html
     #   0-----3-----4-----0
     #   | \   | \   | \   | \   |
-    #   |   \ |   \ |    \|   \ | 
+    #   |   \ |   \ |    \|   \ |
     #   1-----8-----7-----1
     #   | \   | \   | \   | \   |
     #   |   \ |   \ |   \ |   \ |
@@ -418,50 +424,52 @@ def test_persistence_intervals_in_dimension():
     #   |   \ |   \ |   \ |   \ |
     #   0-----3-----4-----0
     st = SimplexTree()
-    st.insert([0,1,8])
-    st.insert([0,3,8])
-    st.insert([3,7,8])
-    st.insert([3,4,7])
-    st.insert([1,4,7])
-    st.insert([0,1,4])
-    st.insert([1,2,5])
-    st.insert([1,5,8])
-    st.insert([5,6,8])
-    st.insert([6,7,8])
-    st.insert([2,6,7])
-    st.insert([1,2,7])
-    st.insert([0,2,3])
-    st.insert([2,3,5])
-    st.insert([3,4,5])
-    st.insert([4,5,6])
-    st.insert([0,4,6])
-    st.insert([0,2,6])
+    st.insert([0, 1, 8])
+    st.insert([0, 3, 8])
+    st.insert([3, 7, 8])
+    st.insert([3, 4, 7])
+    st.insert([1, 4, 7])
+    st.insert([0, 1, 4])
+    st.insert([1, 2, 5])
+    st.insert([1, 5, 8])
+    st.insert([5, 6, 8])
+    st.insert([6, 7, 8])
+    st.insert([2, 6, 7])
+    st.insert([1, 2, 7])
+    st.insert([0, 2, 3])
+    st.insert([2, 3, 5])
+    st.insert([3, 4, 5])
+    st.insert([4, 5, 6])
+    st.insert([0, 4, 6])
+    st.insert([0, 2, 6])
     st.compute_persistence(persistence_dim_max=True)
-    
+
     H0 = st.persistence_intervals_in_dimension(0)
-    assert np.array_equal(H0, np.array([[ 0., float("inf")]]))
+    assert np.array_equal(H0, np.array([[0.0, float("inf")]]))
     H1 = st.persistence_intervals_in_dimension(1)
-    assert np.array_equal(H1, np.array([[ 0., float("inf")], [ 0., float("inf")]]))
+    assert np.array_equal(H1, np.array([[0.0, float("inf")], [0.0, float("inf")]]))
     H2 = st.persistence_intervals_in_dimension(2)
-    assert np.array_equal(H2, np.array([[ 0., float("inf")]]))
+    assert np.array_equal(H2, np.array([[0.0, float("inf")]]))
     # Test empty case
     assert st.persistence_intervals_in_dimension(3).shape == (0, 2)
 
+
 def test_equality_operator():
     st1 = SimplexTree()
     st2 = SimplexTree()
 
     assert st1 == st2
 
-    st1.insert([1,2,3], 4.)
+    st1.insert([1, 2, 3], 4.0)
     assert st1 != st2
 
-    st2.insert([1,2,3], 4.)
+    st2.insert([1, 2, 3], 4.0)
     assert st1 == st2
 
+
 def test_simplex_tree_deep_copy():
     st = SimplexTree()
-    st.insert([1, 2, 3], 0.)
+    st.insert([1, 2, 3], 0.0)
     # compute persistence only on the original
     st.compute_persistence()
 
@@ -480,14 +488,15 @@ def test_simplex_tree_deep_copy():
 
     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_deep_copy_constructor():
     st = SimplexTree()
-    st.insert([1, 2, 3], 0.)
+    st.insert([1, 2, 3], 0.0)
     # compute persistence only on the original
     st.compute_persistence()
 
@@ -506,56 +515,132 @@ def test_simplex_tree_deep_copy_constructor():
 
     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_constructor_exception():
     with pytest.raises(TypeError):
-        st = SimplexTree(other = "Construction from a string shall raise an exception")
+        st = SimplexTree(other="Construction from a string shall raise an exception")
+
+
+def test_create_from_array():
+    a = np.array([[1, 4, 13, 6], [4, 3, 11, 5], [13, 11, 10, 12], [6, 5, 12, 2]])
+    st = SimplexTree.create_from_array(a, max_filtration=5.0)
+    assert list(st.get_filtration()) == [([0], 1.0), ([3], 2.0), ([1], 3.0), ([0, 1], 4.0), ([1, 3], 5.0)]
+
+
+def test_insert_edges_from_coo_matrix():
+    try:
+        from scipy.sparse import coo_matrix
+        from scipy.spatial import cKDTree
+    except ImportError:
+        print("Skipping, no SciPy")
+        return
+
+    st = SimplexTree()
+    st.insert([1, 2, 7], 7)
+    row = np.array([2, 5, 3])
+    col = np.array([1, 4, 6])
+    dat = np.array([1, 2, 3])
+    edges = coo_matrix((dat, (row, col)))
+    st.insert_edges_from_coo_matrix(edges)
+    assert list(st.get_filtration()) == [
+        ([1], 1.0),
+        ([2], 1.0),
+        ([1, 2], 1.0),
+        ([4], 2.0),
+        ([5], 2.0),
+        ([4, 5], 2.0),
+        ([3], 3.0),
+        ([6], 3.0),
+        ([3, 6], 3.0),
+        ([7], 7.0),
+        ([1, 7], 7.0),
+        ([2, 7], 7.0),
+        ([1, 2, 7], 7.0),
+    ]
+
+    pts = np.random.rand(100, 2)
+    tree = cKDTree(pts)
+    edges = tree.sparse_distance_matrix(tree, max_distance=0.15, output_type="coo_matrix")
+    st = SimplexTree()
+    st.insert_edges_from_coo_matrix(edges)
+    assert 100 < st.num_simplices() < 1000
+
+
+def test_insert_batch():
+    st = SimplexTree()
+    # vertices
+    st.insert_batch(np.array([[6, 1, 5]]), np.array([-5.0, 2.0, -3.0]))
+    # triangles
+    st.insert_batch(np.array([[2, 10], [5, 0], [6, 11]]), np.array([4.0, 0.0]))
+    # edges
+    st.insert_batch(np.array([[1, 5], [2, 5]]), np.array([1.0, 3.0]))
+
+    assert list(st.get_filtration()) == [
+        ([6], -5.0),
+        ([5], -3.0),
+        ([0], 0.0),
+        ([10], 0.0),
+        ([0, 10], 0.0),
+        ([11], 0.0),
+        ([0, 11], 0.0),
+        ([10, 11], 0.0),
+        ([0, 10, 11], 0.0),
+        ([1], 1.0),
+        ([2], 1.0),
+        ([1, 2], 1.0),
+        ([2, 5], 4.0),
+        ([2, 6], 4.0),
+        ([5, 6], 4.0),
+        ([2, 5, 6], 4.0),
+    ]
+
 
 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 = 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)
 
     def blocker(simplex):
         try:
             # Block all simplices that contain vertex 6
             simplex.index(6)
-            print(simplex, ' is blocked')
+            print(simplex, " is blocked")
             return True
         except ValueError:
-            print(simplex, ' is accepted')
-            st.assign_filtration(simplex, st.filtration(simplex) + 1.)
+            print(simplex, " is accepted")
+            st.assign_filtration(simplex, st.filtration(simplex) + 1.0)
             return False
 
     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]) == 4.
-    assert st.filtration([0,1,3]) == 5.
-    assert st.filtration([0,2,3]) == 6.
-    assert st.filtration([1,2,3]) == 6.
+    assert st.find([4, 5, 6]) == False
+    assert st.filtration([0, 1, 2]) == 4.0
+    assert st.filtration([0, 1, 3]) == 5.0
+    assert st.filtration([0, 2, 3]) == 6.0
+    assert st.filtration([1, 2, 3]) == 6.0
 
     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]) == 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]) == 7.
+    assert st.find([4, 5, 6]) == False
+    assert st.filtration([0, 1, 2]) == 4.0
+    assert st.filtration([0, 1, 3]) == 5.0
+    assert st.filtration([0, 2, 3]) == 6.0
+    assert st.filtration([1, 2, 3]) == 6.0
+    assert st.filtration([0, 1, 2, 3]) == 7.0