Always treat bars with infinite death as essential bars

[ulupo]

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 #29 (comment)

[ulupo]

A preliminary approach is in ulupo#1

[MonkeyBreaker]

Maybe we could add the diamond example and the output when running in sparse and dense, where we can observe an unexpected difference ?

[ulupo]

@MonkeyBreaker indeed. Here goes:

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

[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

(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

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

(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)])

[MonkeyBreaker]

Perfect, I will just edit your comment to add the python language highlight