[MRG] Fix segfault in ot.emd when no cost is finite

RELEASES.md

@@ -15,6 +15,7 @@
 
 #### Closed issues
 
+- Fix `ot.emd`, `ot.emd2` and the solvers built on them (e.g. `ot.solve`, Gromov-Wasserstein) crashing the interpreter with a segmentation fault when no entry of the cost matrix is finite; the problem is now reported as infeasible (PR #869)
 - Remove a leftover debug `print` from `ot.utils.projection_sparse_simplex` with `axis=1`, and make the `ot.datasets.make_gauss_hd` docstring a raw string so importing `ot` no longer emits a `SyntaxWarning` (PR #860)
 - Fix `ot.dist` ignoring the weights `w` for `metric="cityblock"`, which returned the unweighted distance although the weights are documented for this metric (PR #859)
 - Fix swapped arguments to `div_to_product` in `ot.gromov.fused_unbalanced_across_spaces_cost`: with `reg_type="independent"` (UCOOT) the entropic terms used the plan marginals as the reference measures and vice versa (PR #855, Issue #854)

ot/lp/_network_simplex.py

@@ -218,6 +218,20 @@ def _prepare_warmstart(potentials_init, need_filter, row_mask, col_mask):
     return alpha_init, beta_init
 
 
+def _emd_c_dense(a, b, M, numItermax, numThreads, alpha_init, beta_init):
+    """Call :any:`emd_c`, except for a cost matrix with no finite entry
+
+    The network simplex crashes the interpreter when every cost is NaN or
+    infinite. Such a problem is reported as infeasible instead, which is what
+    the solver already returns when only some of the costs are not finite.
+    """
+    if M.size > 0 and not np.isfinite(M).any():
+        n1, n2 = M.shape
+        infeasible = 0  # ProblemType.INFEASIBLE in emd_wrap
+        return np.zeros((n1, n2)), 0.0, np.zeros(n1), np.zeros(n2), infeasible
+    return emd_c(a, b, M, numItermax, numThreads, alpha_init, beta_init)
+
+
 def emd(
     a,
     b,
@@ -459,7 +473,7 @@ def emd(
         )
 
         # Dense solver
-        G, cost, u, v, result_code = emd_c(
+        G, cost, u, v, result_code = _emd_c_dense(
             a_solver, b_solver, M_solver, numItermax, numThreads, alpha_init, beta_init
         )
 
@@ -776,7 +790,7 @@ def f(b):
             )
 
             # Solve dense EMD
-            G, cost, u, v, result_code = emd_c(
+            G, cost, u, v, result_code = _emd_c_dense(
                 a_solver,
                 b_solver,
                 M_solver,

test/test_ot.py

@@ -266,6 +266,28 @@ def test_emd_empty():
     np.testing.assert_allclose(w, 0)
 
 
+@pytest.mark.parametrize("value", [np.nan, np.inf])
+def test_emd_no_finite_cost(value):
+    # the network simplex crashed the interpreter when no cost was finite;
+    # the problem is reported as infeasible, as for partly non-finite costs
+    a = ot.utils.unif(4)
+    b = ot.utils.unif(3)
+    M = np.full((4, 3), value)
+
+    with pytest.warns(UserWarning, match="infeasible"):
+        G, log = ot.emd(a, b, M, log=True)
+    np.testing.assert_allclose(G, np.zeros((4, 3)))
+    assert log["result_code"] == 0
+
+    with pytest.warns(UserWarning, match="infeasible"):
+        ot.emd2(a, b, M)
+
+    # Gromov-Wasserstein builds such a cost matrix from a NaN structure
+    C1 = np.full((4, 4), value)
+    C2 = np.ones((3, 3))
+    assert np.isnan(ot.gromov.gromov_wasserstein2(C1, C2, a, b))
+
+
 def test_emd2_multi():
     n = 500  # nb bins