Failed assertion when calling dionysus.wasserstein_distance for largish diagrams

[thsutton]

I'm attempting to use dionysus.wasserstein_distance to compare a bunch of diagrams but have run into a crash due to a failed assertion in the C++ code.

My Python code looks a bit like this (after I remove a bunch of control flow, etc):

filtration1 = dionysus.fill_rips(dist1, K, R)
filtration2 = dionysus.fill_rips(dist2, K, R)
homology1 = dionysus.homology_persistence(filtration1)
homology2 = dionysus.homology_persistence(filtration2)
data[i] = dionysus.init_diagrams(homology1, filtration1)
data[j] = dionysus.init_diagrams(homology2, filtration2)
x = data[i][dim]
y = data[j][dim]
print("k=%d, |V[%d]|=%d, |V[%d]|=%d" % (dim, i, len(x), j, len(y)))
d = dionysus.wasserstein_distance(x, y)
logging.debug("\t\t\t\td(%d, %d) = %s", i, j, d)

....
k=2, |V[1]|=2300, |V[0]|=5984
                d(1, 0) = inf
k=2, |V[2]|=1540, |V[0]|=5984
                d(2, 0) = inf
k=2, |V[2]|=1540, |V[1]|=2300
                d(2, 1) = inf
k=2, |V[3]|=5984, |V[0]|=5984
Assertion failed: (dnn_points_.size() < _items.size()), function AuctionOracleKDTreeRestricted, file /Users/thsutton/Documents/TDA/assignment/dionysus/ext/hera/wasserstein/include/auction_oracle_kdtree_restricted.hpp, line 114.

This is an assignment but I can make the code and data available by privately if required.

[thsutton]

Oh, and it's a debug build from 1688d12c0d327a3bdd31489df9abbcd479384067

[anigmetov]

Hello, for now it seems that the problem is in my code, in Hera. Thanks for reporting. Can you send the diagrams x and y that cause the failure to me ( nigmetov e-mail symbol tugraz dot at)? You can do #####################################

f_x = open("dgm_x.txt", "w")
# I don't know exactly how to extract coords from x, y, something like that
# might work. If not, I hope it should be easy to modify the next line appropriately
for a, b in x:
    f_x.write("%f %f\n" % (a, b))
f_x.close()
f_y = open("dgm_y.txt", "w")
for a, b in y:
    f_y.write("%f %f\n" % (a, b))
f_y.close()

##################################### to save the data to the files dgm_x.txt and dgm_y.txt, and I only need these two files.

[thsutton]

I emailed files of the two failing diagrams to @grey-narn.

[anigmetov]

OK, the fix was easy: the diagrams have the same number of points at infinity and no finite points, so I added if (A.empty() and B.empty()) return 0.0; to wasserstein_cost_vec function before calling auction, and pushed the commit to the public Hera repo (added several test cases, too). Dmitry, now you can pull the code to Dionysus. Thanks again to @thsutton for reporting.

[mrzv]

Thanks, @grey-narn!