Improved tests and docs for wasserstein_1d

Author: rtavenar

ot/__init__.py

@@ -22,7 +22,7 @@
 from . import stochastic
 
 # OT functions
-from .lp import emd, emd2, emd_1d, emd2_1d
+from .lp import emd, emd2, emd_1d, emd2_1d, wasserstein_1d, wasserstein2_1d
 from .bregman import sinkhorn, sinkhorn2, barycenter
 from .da import sinkhorn_lpl1_mm
 
@@ -32,5 +32,6 @@
 __version__ = "0.5.1"
 
 __all__ = ["emd", "emd2", 'emd_1d', "sinkhorn", "sinkhorn2", "utils", 'datasets',
-           'bregman', 'lp', 'tic', 'toc', 'toq', 'gromov', 'emd_1d', 'emd2_1d',
+           'bregman', 'lp', 'tic', 'toc', 'toq', 'gromov',
+           'emd_1d', 'emd2_1d', 'wasserstein_1d', 'wasserstein2_1d',
            'dist', 'unif', 'barycenter', 'sinkhorn_lpl1_mm', 'da', 'optim']

ot/lp/__init__.py

@@ -530,13 +530,13 @@ def emd2_1d(x_a, x_b, a=None, b=None, metric='sqeuclidean', p=1., dense=True,
 
 
 def wasserstein_1d(x_a, x_b, a=None, b=None, p=1., dense=True, log=False):
-    """Solves the Wasserstein distance problem between 1d measures and returns
+    """Solves the p-Wasserstein distance problem between 1d measures and returns
     the OT matrix
 
 
     .. math::
-        \gamma = arg\min_\gamma \left(\sum_i \sum_j \gamma_{ij}
-            |x_a[i] - x_b[j]|^p \right)^{1/p}
+        \gamma = arg\min_\gamma \left( \sum_i \sum_j \gamma_{ij}
+            |x_a[i] - x_b[j]|^p \\right)^{1/p}
 
         s.t. \gamma 1 = a,
              \gamma^T 1= b,
@@ -617,15 +617,14 @@ def wasserstein_1d(x_a, x_b, a=None, b=None, p=1., dense=True, log=False):
                   dense=dense, log=log)
 
 
-def wasserstein2_1d(x_a, x_b, a=None, b=None, metric='sqeuclidean', p=1.,
-                    dense=True, log=False):
-    """Solves the Wasserstein distance problem between 1d measures and returns
+def wasserstein2_1d(x_a, x_b, a=None, b=None, p=1., dense=True, log=False):
+    """Solves the p-Wasserstein distance problem between 1d measures and returns
     the loss
 
 
     .. math::
         \gamma = arg\min_\gamma \left( \sum_i \sum_j \gamma_{ij}
-            |x_a[i] - x_b[j]|^p \right)^{1/p}
+            |x_a[i] - x_b[j]|^p \\right)^{1/p}
 
         s.t. \gamma 1 = a,
              \gamma^T 1= b,

ot/lp/emd_wrap.pyx

@@ -13,6 +13,7 @@ cimport numpy as np
 from ..utils import dist
 
 cimport cython
+cimport libc.math as math
 
 import warnings
 
@@ -159,7 +160,7 @@ def emd_1d_sorted(np.ndarray[double, ndim=1, mode="c"] u_weights,
         elif metric == 'cityblock' or metric == 'euclidean':
             m_ij = abs(u[i] - v[j])
         elif metric == 'minkowski':
-            m_ij = abs(u[i] - v[j]) ** p
+            m_ij = math.pow(abs(u[i] - v[j]), p)
         else:
             m_ij = dist(u[i].reshape((1, 1)), v[j].reshape((1, 1)),
                         metric=metric)[0, 0]

test/test_ot.py

@@ -85,6 +85,29 @@ def test_emd_1d_emd2_1d():
     np.testing.assert_raises(AssertionError, ot.emd_1d, u, v, [], [])
 
 
+def test_wass_1d():
+    # test emd1d gives similar results as emd
+    n = 20
+    m = 30
+    rng = np.random.RandomState(0)
+    u = rng.randn(n, 1)
+    v = rng.randn(m, 1)
+
+    M = ot.dist(u, v, metric='sqeuclidean')
+
+    G, log = ot.emd([], [], M, log=True)
+    wass = log["cost"]
+
+    G_1d, log = ot.wasserstein_1d(u, v, [], [], p=2., log=True)
+    wass1d = log["cost"]
+
+    # check loss is similar
+    np.testing.assert_allclose(np.sqrt(wass), wass1d)
+
+    # check G is similar
+    np.testing.assert_allclose(G, G_1d)
+
+
 def test_emd_empty():
     # test emd and emd2 for simple identity
     n = 100