New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Proposal for dealing with explicit infinite deaths #1

Closed

ulupo wants to merge 55 commits into feature/return_flag_persistence_generators from feature/return_flag_persistence_generators_temp

Owner

ulupo commented Oct 31, 2021 •

edited

Even in the case of dense matrices, any bar which has infinity as death must be considered an essential bar. Similarly, for the dimension-0 computation, we must treat an edge with value infinity as non-existent.

We did not notice this so far because we were only printing barcodes, but this mismatch gives mathematically wrong results when we want to return generators.

This is a proposal that does not yet yield the correct results in higher dimensions (and I don't know why).

For example,

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([[0, 0, 0, 0]], dtype=int64)],
 array([0], dtype=int64),
 [array([], shape=(0, 2), dtype=int64), array([], shape=(0, 2), dtype=int64)])

so I failed to create an essential generator in dimension 2 and created a corrupt finite one instead.

ulupo and others added 30 commits

September 24, 2021 15:52


          WIP first sketch (won't compile)

a871bcb


          [CPP] Fix compilation issue, refactor retrieving representative simpl…

09ddc8d

…ices

Currently it should not work

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Retrieve representative simplices from computation in CPP

ec5ca78

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Disable sorting of barcodes to keep alignment with representati…

…ve simplices

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Fix not enough memory pre-allocated for all dimensions

32f23fd

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [PY] Add draft of Python interface and Bindings support

2efdff6

for representative simplices

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Retrieve only "True" representative simplices

ae208e6

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Retrieve representative simplices with an indexing 1:1 with bar…

15bc664

…codes

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Refactor with clang-format

922d58b

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Update size of representative simplices, to go from dim 1 to dim N

e4f302a

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Allow to sort barcodes and reproduce the same sort on rep simpl…

eb8426d

…ices

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Add representative simplicies for dim 0

d244359

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [PY] Add output of representative simplices for dim 0 (finite and ess)

567d75a

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [TEST] Add test for verify that all edges in rpsm are valid in the dm

9e963c1

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          Remove array conversions, enforce int type for indices, ensure array …

9c3bc46

…shapes are correct


          [PY] Change dtype to int64 in flag generators

780c0fc


          [CPP] Refactor get_youngest_edge_simplex function

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Refactor code related to representative simplicies

1e69e96

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Fix vertices order in dimension 0 for representative simplicies

fc087b5

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [PY/TEST] Add exception if user uses collapser + representatives.

e7b71cf

test added to verify this behavior

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [DOC] Add draft for input and output documentation of representative

fbe1209

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [EXAMPLE] Add a section about representative simplices example

c7f5d46

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          Fix birth simplices for finite 0-dimensional bars

7a81203


          Actual fix for 0-dimensional persistence generator

cbb1c83


          [CPP] Refactor @ulupo fix for representative in dim 0

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          Add @ctralie in copyright notice

2d58ce9


          Add inline comment as suggested by @MonkeyBreaker

6deb759


          [CPP] Update as suggested by @ulupo method link of class union_find

7a57c76

The link function now also returns the birth vertex and allows to simplify logic in the computation of 0-dim

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Revert latest commit

b085aaf

link method works on the vertex index and the rank of the vertices
but the logic to retrieve the correct birth (value_t) value requires comparison of birth values
that IMO should not be integrated in the link method. To be discussed later
But for the moment I prefer revert this change

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Attempt fix for link_and_get_birth_vertex

c7bab16

- Extract birth using dset.get_birth
- FIXME avoid calling find twice for u and v by using an "unsafe" version of link_and_get_birth_vertex called _link_and_get_birth_vertex

ulupo and others added 15 commits

October 26, 2021 17:13


          Update kwarg name in tests

d9fb0ba


          Change name to return_generators in bindings

67b3687


          Mathematically intelligible solution to 0-dimensional generators comp…

1b02363

…utation


          Improve jupyter tutorial

a529ea3


          Use np.inf instead of np.infty

5bb87a3


          Notebook refactor

5d12970


          Mention of ripser.py in basic tutorial

3c355f0


          [CPP] Add inline comments in link_roots_and_get_birth_vertex

61a3e69


          [TEST] add to test_..._diagonal_dim0 check for the number of points

9a40889

Fix that no value in the diagonal should be bigger than a non diagonal value

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [TEST] Add simple test for finite generate at dim > 0 suggested by @u…

0ea9b1e

…lupo

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Rename representative by generator

fe4c0fb

Refactor slighly compute_dim_0

Signed-off-by: julian <julian.burellaperez@heig-vd.ch>


          [CPP] Fix typos in two inline comments

bbaacd5


          Merge branch 'giotto-ai:main' into feature/return_flag_persistence_ge…

67664dd

…nerators


          [CPP] fix small errors in inline comments

dc01c53


          Proposal for dealing with explicit infinite deaths

c3520f9

ulupo requested a review from MonkeyBreaker

October 31, 2021 10:05

ulupo added 2 commits

October 31, 2021 14:57


          Fix previous attempt, seems to be working now

777a788


          Add inline comment

371fd71

ulupo commented

View reviewed changes

gph/src/ripser.h

@@ @@ -1562,11 +1585,9 @@ class ripser @@
                                     std::greater<>());
               #endif
                           for (size_t i = 0; i < idx_essential; ++i) {
-                              if (!std::isinf(essential_pair[i])) {

Owner Author

ulupo Oct 31, 2021

This check is redundant now that I explicitly compute is_born

ulupo commented

View reviewed changes

gph/src/ripser.h Outdated

Comment on lines 1488 to 1490

+                                  // TODO: these will need special attention, if output
+                                  // happens after the reduction, not during
+                                  break;

Owner Author

ulupo Oct 31, 2021

I never really understood this comment but I wonder if it still makes sense in this new position

ulupo commented

View reviewed changes

gph/src/ripser.h

@@ @@ -1340,6 +1342,7 @@ class ripser @@
                           }
                           diameter_entry_t e;
+                          bool is_essential = false;

Owner Author

ulupo Oct 31, 2021

Performance impact? (I'd say insignificant?)

ulupo commented

View reviewed changes

gph/src/ripser.h

-                                          flag_persistence_generators.finite_0.push_back(
-                                              {birth_vertex.second, std::max(v1, v2),
-                                               std::min(v1, v2)});
+                          if (get_diameter(e) != std::numeric_limits<value_t>::infinity()) {

Owner Author

ulupo Oct 31, 2021

Again, if an edge is infinite it means it's not there, so we should not use it to link vertices.

ulupo commented

View reviewed changes

gph/src/ripser.h

Comment on lines +1443 to +1450

+                                          /* Pairs should only be extracted if insertion was
+                                           * first one! */
+                                          size_t location = last_diameter_index++;
+                                          diameters[location] = death;
+                                          auto first_ins =
+                                              deaths
+                                                  .insert({get_index(get_entry(pivot)), location})
+                                                  .second;

Owner Author

ulupo Oct 31, 2021

I don't yet understand these lines too well so please check that them now being inside this if cannot lead to wrong results.

Collaborator

MonkeyBreaker Oct 31, 2021

It's an explanation of what is going on. Because the reduction is done in parallel. You can have at a certain point 2 differents workers getting the same pivot as the last reduction. If that happens, we ensure with this check that only one barcode is produced instead of two.

I think that the comment should remain for anyone looking at the code

Owner Author

ulupo Oct 31, 2021

Sure, I wasn't complaining about the comment, just wondering about correctness of the code in this part especially.

Collaborator

MonkeyBreaker Oct 31, 2021

Oh !
Is it more clear with my previous explanation ?

ulupo commented

View reviewed changes

gph/src/ripser.h

Comment on lines +1455 to +1456

		if (first_ins) {
		if (return_flag_persistence_generators) {

Owner Author

ulupo Oct 31, 2021

Could be dealt with with an && as it was before

ulupo commented

View reviewed changes

Owner Author

ulupo left a comment

Global comment: technically we should be checking for positive infinity only, I'm not sure whether this is what's happening now with the == syntax

ulupo mentioned this pull request

Always treat bars with infinite death as essential bars giotto-ai/giotto-ph#37

Closed

MonkeyBreaker force-pushed the feature/return_flag_persistence_generators branch from 25fd98f to 0abb23b Compare

October 31, 2021 18:23

ulupo deleted the branch feature/return_flag_persistence_generators

November 1, 2021 10:07

ulupo closed this

ulupo deleted the feature/return_flag_persistence_generators_temp branch

November 26, 2021 17:17

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment