Erroneous assertion?

[gspr]

This assertion seems like it can fail in a lot of normal (non-error) circumstances. Perhaps I'm reading the code wrong – and I apologize for the noise if so – but won't this strict inequality fail any time a persistence diagram contains only points satisfying the is_normal test in the loop above?

[anigmetov]

It will only fail if the other diagram is empty: for every normal, i.e., off-diagonal, point I add a diagonal projection of this point to the other diagram. User input is not supposed to have those. I think I have some shortcut that checks for empty diagram case before I actually get to the auction (in that case everything is matched to the diagonal and the distance is (the root of) the sum of persistences of every point of the non-empty diagram raised to the appropriate power).

[gspr]

It will only fail if the other diagram is empty: for every normal, i.e., off-diagonal, point I add a diagonal projection of this point to the other diagram. User input is not supposed to have those. I think I have some shortcut that checks for empty diagram case before I actually get to the auction (in that case everything is matched to the diagonal and the distance is (the root of) the sum of persistences of every point of the non-empty diagram raised to the appropriate power).

Thanks for the super quick reply! Maybe I'm confused, but this does seem to crop up in a test in GUDHI.

[anigmetov]

I looked at the failing tests; the problem was that I did not consider diagrams that only include points at infinity as empty. I modified the logic accordingly, so the execution now should not reach the assertion for diag9, diag10. Can you pull the latest commit and verify that your tests pass now?

[gspr]

I'm still seeing the same GUDHI test fail. Are you seeing it pass?

[anigmetov]

We are talking about diag9, diag10, aren't we? I modified my reader so that the points with the same x and y coordinate are not ignored, and no assertion fires (and if I put assert(false), somewhere, this one gets triggered, so assertions are working):

hera/build/wasserstein git:(master) x cat a.txt 
-inf -inf
inf inf
hera/build/wasserstein git:(master) x cat b.txt 
-inf -inf
inf inf
0 1
hera/build/wasserstein git:(master) x ./wasserstein_dist a.txt b.txt 
0.5

[gspr]

On 6 December 2024 01:27:56 CET, Arnur Nigmetov ***@***.***> wrote: We are talking about diag9, diag10, aren't we? I modified my reader so that the points with the same x and y coordinate are not ignored, and no assertion fires (and if I put assert(false), somewhere, this one gets triggered, so assertions are working): ``` hera/build/wasserstein git:(master) x cat a.txt -inf -inf inf inf hera/build/wasserstein git:(master) x cat b.txt -inf -inf inf inf 0 1 hera/build/wasserstein git:(master) x ./wasserstein_dist a.txt b.txt 0.5 ```

My apologies - the commit does fix the assertion failure that I reported in this issue. However, that same test in GUDHI now triggers a different assertion failure, which confused me. I haven't had time to figure out where that is, or whether it's on the Hera or the GUDHI side. Thus: thanks for the fix, and feel free to close the issue. I'll open a new one if the new fault is actually in Hera. -- Sent from my Android device with K-9 Mail. Please excuse my brevity.

[VincentRouvreau]

I confirm your commit fixed the assertion.

The assertion that is raised now comes from GUDHI : https://github.com/GUDHI/gudhi-devel/blob/350a3dbd775e811f75b19903fd319978a4f5c197/src/python/gudhi/hera/wasserstein.cc#L83
when doing:

from gudhi.hera import wasserstein_distance
import numpy as np
diag9 = np.array([[-np.inf, -np.inf], [np.inf, np.inf]])
diag10 = np.array([[0,1], [-np.inf, -np.inf], [np.inf, np.inf]])
wasserstein_distance(diag9, diag10, matching=True, internal_p=1., order=1.)
# RuntimeError: Caught an unknown exception!
# n1=2 - n2=3

For me this is a GUDHI bug as we set

  int n1 = boost::size(diag1);
  int n2 = boost::size(diag2);
  ...
  GUDHI_CHECK(n1==n2, "unexpected bug in Hera?");

And the next GUDHI_CHECK(v1[i][0]==v2[i][0] && v1[i][1]==v2[i][1], "unexpected bug in Hera?"); in the loop is also failing. Maybe this is because the comment // bug in Hera, matching_a_to_b_ is empty if on(diag9, diag10, matching=True, internal_p=1., order=1.)e diagram is empty or both diagrams contain the same points is no more true. @mglisse ?

If I remove the if(res.matching_a_to_b_.size() == 0) { ... }, part, it returns:

wasserstein_distance(diag9, diag10, matching=True, internal_p=1., order=1.)
(1.0, array([[-1,  1],
       [ 1, -1],
       [-1,  2],
       [ 0, -1]], dtype=int32))

Which is not exactly what I have with POT version:

wasserstein_distance_pot(diag9, diag10, matching=True, internal_p=1., order=1.)
(1.0, array([[-1,  0],
       [ 0, -1],
       [ 1, -1],
       [-1,  1],
       [-1,  2]]))

[mglisse]

It looks like in this part of the gudhi workaround, I fail to handle points on the diagonal properly (they are handled later), I'll have to fix that.
Note that this code in gudhi is needed because Hera returns an empty matching_a_to_b_ for non-empty diagrams...

[mglisse]

hera/build/wasserstein git:(master) x cat a.txt 
-inf -inf
inf inf
hera/build/wasserstein git:(master) x cat b.txt 
-inf -inf
inf inf
0 1
hera/build/wasserstein git:(master) x ./wasserstein_dist a.txt b.txt 
0.5

If I ask for the matching, I get something empty. I would expect at the very least the point "0 1" to be matched to the diagonal, but the logic for empty finite diagrams seems to ignore the matching. It could also make sense to make sure that every single input point is matched to something (although I can understand if you decide to skip the diagonal points).

#include <iostream>
#include <limits>
#include <hera/wasserstein.h>
const double inf = std::numeric_limits<double>::infinity();
typedef hera::DiagramPoint<double> P; // { x, y, idx }
int main(){
  std::vector<P> diag1 {{-inf,-inf,0},{inf,inf,1}};
  std::vector<P> diag2 {{-inf,-inf,0},{inf,inf,1},{0,1,2}};
  hera::AuctionParams<double> params;
  params.return_matching = true;
  params.match_inf_points = true;
  auto res = hera::wasserstein_cost_detailed(diag1, diag2, params);
  std::cout << "matching_a_to_b_\n";
  for (auto&p : res.matching_a_to_b_)
    std::cout << p.first << '\t' << p.second << '\n';
  std::cout << "matching_b_to_a_\n";
  for (auto&p : res.matching_b_to_a_)
    std::cout << p.first << '\t' << p.second << '\n';
}

matching_a_to_b_
matching_b_to_a_

[mglisse]

I also get an empty matching for 2 identical diagrams

  std::vector<P> diag1 {{0,2,0},{0,inf,1},{1,3,2}};
  std::vector<P> diag2 {{0,2,0},{0,inf,1},{1,3,2}};

(I am hoping that instead of having to fix my workarounds in Gudhi, you can directly fix Hera...)