Repository: giotto-ai/giotto-ph (language: python)

Issue:
Always treat bars with infinite death as essential bars
During work on this PR #29, @MonkeyBreaker and I realized that the current core logic in the matrix reduction algorithm makes it so that we break the following promise to the user: that passing a dense matrix with a few `np.inf` will lead to exactly the same output as passing a sparse version where the `np.inf` entries are just missing entries.  This was never noticed at the level of barcodes because a "finite" bar with death equal to infinity looks just the same as an infinite bar with the same birth.  But flag generators don't lie: we now end up telling the user there are never any generators for essential classes if the matrix is dense, no matter what the matrix is or represents.

Because changing this requires potentially delicate changes to the reduction routine, we decided to discuss what to do about this separately, and open a separate PR if necessary after this one is merged.

_Originally posted by @ulupo in https://github.com/giotto-ai/giotto-ph/issues/29#issuecomment-955710360_
Always treat bars with infinite death as essential bars
During work on this PR #29, @MonkeyBreaker and I realized that the current core logic in the matrix reduction algorithm makes it so that we break the following promise to the user: that passing a dense matrix with a few `np.inf` will lead to exactly the same output as passing a sparse version where the `np.inf` entries are just missing entries.  This was never noticed at the level of barcodes because a "finite" bar with death equal to infinity looks just the same as an infinite bar with the same birth.  But flag generators don't lie: we now end up telling the user there are never any generators for essential classes if the matrix is dense, no matter what the matrix is or represents.

Because changing this requires potentially delicate changes to the reduction routine, we decided to discuss what to do about this separately, and open a separate PR if necessary after this one is merged.

_Originally posted by @ulupo in https://github.com/giotto-ai/giotto-ph/issues/29#issuecomment-955710360_


Discussion / hints:
A preliminary approach is in https://github.com/ulupo/giotto-ph/pull/1
Maybe we could add the diamond example and the output when running in sparse and dense, where we can observe an unexpected difference ?
@MonkeyBreaker indeed. Here goes:

```python
import numpy as np
from gph.python import ripser_parallel

diamond_dm = np.array(
    [[0,      1,      np.inf, 1,      1,      1],
     [0,      0,      1,      np.inf, 1,      1],
     [0,      0,      0,      1,      1,      1],
     [0,      0,      0,      0,      1,      1],
     [0,      0,      0,      0,      0,      np.inf],
     [0,      0,      0,      0,      0,      0]],
    dtype=np.float64
)
diamond_dm += diamond_dm.T

res = ripser_parallel(diamond_dm, metric="precomputed", maxdim=2, return_generators=True)
```
yields the correct barcode
```python
[array([[ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0., inf]]),
 array([], shape=(0, 2), dtype=float64),
 array([[ 1., inf]])]
```
but the incorrect generators
```python
(array([[5, 5, 3],
       [3, 5, 2],
       [2, 5, 1],
       [1, 5, 0],
       [4, 4, 3]], dtype=int64),
 [array([], shape=(0, 4), dtype=int64), array([[1, 0, 5, 4]], dtype=int64)],
 array([0], dtype=int64),
 [array([], shape=(0, 2), dtype=int64), array([], shape=(0, 2), dtype=int64)])
```

The sparse version
```python
from scipy.sparse import coo_matrix

diamond_dm = np.where(diamond_dm == np.inf, 0, diamond_dm)
diamond_dm = coo_matrix(diamond_dm)
```
gives the correct
```python
(array([[5, 5, 3],
       [3, 5, 2],
       [2, 5, 1],
       [1, 5, 0],
       [4, 4, 3]], dtype=int64),
 [array([], shape=(0, 4), dtype=int64), array([], shape=(0, 4), dtype=int64)],
 array([0], dtype=int64),
 [array([], shape=(0, 2), dtype=int64), array([[1, 0]], dtype=int64)])
```
Perfect, I will just edit your comment to add the `python` language highlight
A preliminary approach is in https://github.com/ulupo/giotto-ph/pull/1
Maybe we could add the diamond example and the output when running in sparse and dense, where we can observe an unexpected difference ?
@MonkeyBreaker indeed. Here goes:

```python
import numpy as np
from gph.python import ripser_parallel

diamond_dm = np.array(
    [[0,      1,      np.inf, 1,      1,      1],
     [0,      0,      1,      np.inf, 1,      1],
     [0,      0,      0,      1,      1,      1],
     [0,      0,      0,      0,      1,      1],
     [0,      0,      0,      0,      0,      np.inf],
     [0,      0,      0,      0,      0,      0]],
    dtype=np.float64
)
diamond_dm += diamond_dm.T

res = ripser_parallel(diamond_dm, metric="precomputed", maxdim=2, return_generators=True)
```
yields the correct barcode
```python
[array([[ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0.,  1.],
        [ 0., inf]]),
 array([], shape=(0, 2), dtype=float64),
 array([[ 1., inf]])]
```
but the incorrect generators
```python
(array([[5, 5, 3],
       [3, 5, 2],
       [2, 5, 1],
       [1, 5, 0],
       [4, 4, 3]], dtype=int64),
 [array([], shape=(0, 4), dtype=int64), array([[1, 0, 5, 4]], dtype=int64)],
 array([0], dtype=int64),
 [array([], shape=(0, 2), dtype=int64), array([], shape=(0, 2), dtype=int64)])
```

The sparse version
```python
from scipy.sparse import coo_matrix

diamond_dm = np.where(diamond_dm == np.inf, 0, diamond_dm)
diamond_dm = coo_matrix(diamond_dm)
```
gives the correct
```python
(array([[5, 5, 3],
       [3, 5, 2],
       [2, 5, 1],
       [1, 5, 0],
       [4, 4, 3]], dtype=int64),
 [array([], shape=(0, 4), dtype=int64), array([], shape=(0, 4), dtype=int64)],
 array([0], dtype=int64),
 [array([], shape=(0, 2), dtype=int64), array([[1, 0]], dtype=int64)])
```
Perfect, I will just edit your comment to add the `python` language highlight

Fix patch:
diff --git a/gph/src/ripser.h b/gph/src/ripser.h
index 00d2581c..ec2b1eb5 100644
--- a/gph/src/ripser.h
+++ b/gph/src/ripser.h
@@ -94,6 +94,7 @@ typedef uint16_t coefficient_t;
 using barcodes_t = std::vector<value_t>;
 
 const size_t num_coefficient_bits = 8;
+constexpr value_t inf_value = std::numeric_limits<value_t>::infinity();
 
 // 1L on windows is ALWAYS 32 bits, when on unix systems is pointer size
 static const index_t max_simplex_index =
@@ -394,7 +395,7 @@ struct sparse_distance_matrix {
                              index_diameter_t{j, vertex_births[j]});
         return (neighbor != neighbors[i].end() && get_index(*neighbor) == j)
                    ? get_diameter(*neighbor)
-                   : std::numeric_limits<value_t>::infinity();
+                   : inf_value;
     }
 
     size_t size() const { return neighbors.size(); }
@@ -713,7 +714,7 @@ class ripser
     value_t compute_diameter(const index_t index, const index_t dim) const
     {
         std::vector<index_t> vertices;
-        value_t diam = -std::numeric_limits<value_t>::infinity();
+        value_t diam = -inf_value;
 
         vertices.resize(dim + 1);
         get_simplex_vertices(index, dim, dist.size(), vertices.rbegin());
@@ -980,6 +981,9 @@ class ripser
         value_t birth, death;
         diameter_index_t birth_vertex;
         for (auto e : edges) {
+            if (get_diameter(e) == inf_value)
+                continue;
+
             get_simplex_vertices(get_index(e), 1, n, vertices_of_edge.rbegin());
             index_t v1 = vertices_of_edge[0];
             index_t v2 = vertices_of_edge[1];
@@ -1014,8 +1018,7 @@ class ripser
         for (index_t i = 0; i < n; ++i) {
             if (dset.find(i) == i) {
                 births_and_deaths_by_dim[0].push_back(dset.get_birth(i));
-                births_and_deaths_by_dim[0].push_back(
-                    std::numeric_limits<value_t>::infinity());
+                births_and_deaths_by_dim[0].push_back(inf_value);
                 if (return_flag_persistence_generators) {
                     flag_persistence_generators.essential_0.push_back(
                         dset.get_birth_idx(i));
@@ -1242,7 +1245,7 @@ class ripser
     edge_t
     get_youngest_edge_simplex(const std::vector<index_t>& vertices_simplex)
     {
-        value_t diam = -std::numeric_limits<value_t>::infinity();
+        value_t diam = -inf_value;
         edge_t edge;
 
         for (index_t i = 0; i < vertices_simplex.size(); ++i) {
@@ -1329,7 +1332,9 @@ class ripser
                                false);
                 pivot = get_pivot(working_coboundary);
             }
+
             diameter_entry_t e;
+            bool is_essential = false;
 
             while (true) {
                 if (get_index(pivot) != -1) {
@@ -1423,49 +1428,69 @@ class ripser
                                            // the new pair; could elide an
                                            // extra atomic load
 
-                        /* Pairs should be extracted if insertion was
-                         * first one ! */
-                        const auto new_idx_finite_bar = idx_finite_bar++;
-                        death_diams[new_idx_finite_bar] = get_diameter(pivot);
-                        auto first_ins =
-                            pivot_to_death_idx
-                                .insert({get_index(get_entry(pivot)),
-                                         new_idx_finite_bar})
-                                .second;
-
-                        /* Only insert when it is the first time this bar is
-                         * encountered
-                         */
-                        if (return_flag_persistence_generators && first_ins) {
-                            // Inessential flag persistence generators in dim >
-                            // 0
-                            std::vector<index_t> vertices_birth(dim + 1);
-                            std::vector<index_t> vertices_death(dim + 2);
-
-                            index_t birth_idx =
-                                get_index(get_entry(column_to_reduce));
-                            index_t death_idx = get_index(get_entry(pivot));
-
-                            get_simplex_vertices(birth_idx, dim, n,
-                                                 vertices_birth.rbegin());
-                            get_simplex_vertices(death_idx, dim + 1, n,
-                                                 vertices_death.rbegin());
-
-                            edge_t birth_edge =
-                                get_youngest_edge_simplex(vertices_birth);
-                            edge_t death_edge =
-                                get_youngest_edge_simplex(vertices_death);
-
-                            finite_generator[new_idx_finite_bar] = {
-                                birth_edge.first, birth_edge.second,
-                                death_edge.first, death_edge.second};
+                        // Infinite death indicates an essential bar in disguise
+                        value_t death = get_diameter(pivot);
+                        is_essential = (death == inf_value);
+
+                        if (!is_essential) {
+                            /* Pairs should only be extracted if insertion was
+                             * first one! */
+                            const auto new_idx_finite_bar = idx_finite_bar++;
+                            death_diams[new_idx_finite_bar] = death;
+                            auto first_ins =
+                                pivot_to_death_idx
+                                    .insert({get_index(get_entry(pivot)),
+                                             new_idx_finite_bar})
+                                    .second;
+
+                            /* Only insert when it is the first time this bar is
+                             * encountered
+                             */
+                            if (first_ins &&
+                                return_flag_persistence_generators) {
+                                // Inessential flag persistence generators in
+                                // dim > 0
+                                std::vector<index_t> vertices_birth(dim + 1);
+                                std::vector<index_t> vertices_death(dim + 2);
+
+                                index_t birth_idx =
+                                    get_index(get_entry(column_to_reduce));
+                                index_t death_idx = get_index(get_entry(pivot));
+
+                                get_simplex_vertices(birth_idx, dim, n,
+                                                     vertices_birth.rbegin());
+                                get_simplex_vertices(death_idx, dim + 1, n,
+                                                     vertices_death.rbegin());
+
+                                edge_t birth_edge =
+                                    get_youngest_edge_simplex(vertices_birth);
+                                edge_t death_edge =
+                                    get_youngest_edge_simplex(vertices_death);
+
+                                finite_generator[new_idx_finite_bar] = {
+                                    birth_edge.first, birth_edge.second,
+                                    death_edge.first, death_edge.second};
+                            }
                         }
 
                         break;
                     }
                 } else {
+                    is_essential = true;
+                    break;
+                }
+            }
+
+            if (is_essential) {
+                // TODO: these will need special attention, if output
+                // happens after the reduction, not during
+
+                const value_t birth = get_diameter(column_to_reduce);
+                /* Infinite birth means the bar is never born and should
+                 * not be included. */
+                if (birth != inf_value) {
                     auto idx_ = idx_essential++;
-                    essential_pair[idx_] = get_diameter(column_to_reduce);
+                    essential_pair[idx_] = birth;
 
                     if (return_flag_persistence_generators) {
                         // Essential flag persistence generators in dim > 0
@@ -1483,12 +1508,9 @@ class ripser
                         essential_generator[idx_] = {birth_edge.first,
                                                      birth_edge.second};
                     }
-
-                    // TODO: these will need special attention, if output
-                    // happens after the reduction, not during
-                    break;
                 }
             }
+
             return index_column_to_reduce;
         })
             ;
@@ -1554,11 +1576,8 @@ class ripser
                       std::greater<>());
 #endif
             for (size_t i = 0; i < idx_essential; ++i) {
-                if (!std::isinf(essential_pair[i])) {
-                    births_and_deaths_by_dim[dim].push_back(essential_pair[i]);
-                    births_and_deaths_by_dim[dim].push_back(
-                        std::numeric_limits<value_t>::infinity());
-                }
+                births_and_deaths_by_dim[dim].push_back(essential_pair[i]);
+                births_and_deaths_by_dim[dim].push_back(inf_value);
             }
         }
 


Test patch:
diff --git a/gph/python/test/test_ripser.py b/gph/python/test/test_ripser.py
index ffcc1021..cc87bcc7 100644
--- a/gph/python/test/test_ripser.py
+++ b/gph/python/test/test_ripser.py
@@ -431,3 +431,25 @@ def test_unsupported_coefficient():
     with pytest.raises(ValueError):
         ripser(X, metric='precomputed',
                coeff=gph_ripser.get_max_coefficient_field_supported()+1)
+
+
+def test_infinite_deaths_always_essential():
+    """Regression test for issue #37"""
+    diamond_dm = np.array(
+        [[0,      1,      np.inf, 1,      1,      1],
+         [0,      0,      1,      np.inf, 1,      1],
+         [0,      0,      0,      1,      1,      1],
+         [0,      0,      0,      0,      1,      1],
+         [0,      0,      0,      0,      0,      np.inf],
+         [0,      0,      0,      0,      0,      0]],
+        dtype=np.float64
+    )
+    diamond_dm += diamond_dm.T
+
+    gens = ripser(diamond_dm, metric="precomputed", maxdim=2,
+                  return_generators=True)["gens"]
+
+    gens_fin_dim1 = gens[1][1]
+
+    # With this example no finite generators in dimension 1 shall be found
+    assert len(gens_fin_dim1) == 0