Merge branch 'master' into semirelaxed_gromov

Author: rflamary

RELEASES.md

@@ -7,6 +7,9 @@
 - Added semi-relaxed (Fused) Gromov-Wasserstein solvers + examples (PR #401)
 - Added Generalized Wasserstein Barycenter solver + example (PR #372), fixed graphical details on the example (PR #376)
 - Added Free Support Sinkhorn Barycenter + example (PR #387)
+- New API for OT solver using function `ot.solve` (PR #388)
+- Backend version of `ot.partial` and `ot.smooth`  (PR #388)
+
 
 #### Closed issues
 
@@ -26,6 +29,14 @@ roughly 2^31) (PR #381)
 - Fixed an issue where a pytorch example would throw an error if executed on a GPU (Issue #389, PR #391)
 - Added a work-around for scipy's bug, where you cannot compute the Hamming distance with a "None" weight attribute. (Issue #400, PR #402)
 - Fixed an issue where the doc could not be built due to some changes in matplotlib's API (Issue #403, PR #402)
+- Replaced Numpy C Compiler with Setuptools C Compiler due to deprecation issues (Issue #408, PR #409)
+- Fixed weak optimal transport docstring (Issue #404, PR #410)
+- Fixed error with parameter `log=True`for `SinkhornLpl1Transport` (Issue #412, 
+PR #413)
+- Fixed an issue about `warn` parameter in `sinkhorn2` (PR #417)
+- Fix an issue where the parameter `stopThr` in `empirical_sinkhorn_divergence` was rendered useless by subcalls 
+  that explicitly specified `stopThr=1e-9` (Issue #421, PR #422).
+- Fixed a bug breaking an example where we would try to make an array of arrays of different shapes (Issue #424, PR #425)
 
 
 ## 0.8.2

examples/gromov/plot_gromov_barycenter.py

@@ -110,8 +110,7 @@ def im2mat(img):
             if shapes[nb][i, j] < 0.95:
                 xs[nb].append([j, 8 - i])
 
-xs = np.array([np.array(xs[0]), np.array(xs[1]),
-               np.array(xs[2]), np.array(xs[3])])
+xs = [np.array(xs[s]) for s in range(S)]
 
 ##############################################################################
 # Barycenter computation

ot/__init__.py

@@ -34,6 +34,7 @@
 from . import regpath
 from . import weak
 from . import factored
+from . import solvers
 
 # OT functions
 from .lp import emd, emd2, emd_1d, emd2_1d, wasserstein_1d
@@ -46,7 +47,7 @@
                         gromov_barycenters, fused_gromov_wasserstein, fused_gromov_wasserstein2)
 from .weak import weak_optimal_transport
 from .factored import factored_optimal_transport
-
+from .solvers import solve
 
 # utils functions
 from .utils import dist, unif, tic, toc, toq
@@ -61,5 +62,5 @@
            'sinkhorn_unbalanced2', 'sliced_wasserstein_distance',
            'gromov_wasserstein', 'gromov_wasserstein2', 'gromov_barycenters', 'fused_gromov_wasserstein', 'fused_gromov_wasserstein2',
             'max_sliced_wasserstein_distance', 'weak_optimal_transport',
-            'factored_optimal_transport', 
-           'smooth', 'stochastic', 'unbalanced', 'partial', 'regpath']
+            'factored_optimal_transport', 'solve',
+           'smooth', 'stochastic', 'unbalanced', 'partial', 'regpath', 'solvers']

ot/backend.py

@@ -854,6 +854,21 @@ def sqrtm(self, a):
         """
         raise NotImplementedError()
 
+    def kl_div(self, p, q, eps=1e-16):
+        r"""
+        Computes the Kullback-Leibler divergence.
+
+        This function follows the api from :any:`scipy.stats.entropy`.
+
+        Parameter eps is used to avoid numerical errors and is added in the log.
+
+        .. math::
+             KL(p,q) = \sum_i p(i) \log (\frac{p(i)}{q(i)}+\epsilon)
+
+        See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.entropy.html
+        """
+        raise NotImplementedError()
+
     def isfinite(self, a):
         r"""
         Tests element-wise for finiteness (not infinity and not Not a Number).
@@ -1158,6 +1173,9 @@ def inv(self, a):
     def sqrtm(self, a):
         return scipy.linalg.sqrtm(a)
 
+    def kl_div(self, p, q, eps=1e-16):
+        return np.sum(p * np.log(p / q + eps))
+
     def isfinite(self, a):
         return np.isfinite(a)
 
@@ -1481,6 +1499,9 @@ def sqrtm(self, a):
         L, V = jnp.linalg.eigh(a)
         return (V * jnp.sqrt(L)[None, :]) @ V.T
 
+    def kl_div(self, p, q, eps=1e-16):
+        return jnp.sum(p * jnp.log(p / q + eps))
+
     def isfinite(self, a):
         return jnp.isfinite(a)
 
@@ -1901,6 +1922,9 @@ def sqrtm(self, a):
         L, V = torch.linalg.eigh(a)
         return (V * torch.sqrt(L)[None, :]) @ V.T
 
+    def kl_div(self, p, q, eps=1e-16):
+        return torch.sum(p * torch.log(p / q + eps))
+
     def isfinite(self, a):
         return torch.isfinite(a)
 
@@ -2248,6 +2272,9 @@ def sqrtm(self, a):
         L, V = cp.linalg.eigh(a)
         return (V * self.sqrt(L)[None, :]) @ V.T
 
+    def kl_div(self, p, q, eps=1e-16):
+        return cp.sum(p * cp.log(p / q + eps))
+
     def isfinite(self, a):
         return cp.isfinite(a)
 
@@ -2608,6 +2635,9 @@ def inv(self, a):
     def sqrtm(self, a):
         return tf.linalg.sqrtm(a)
 
+    def kl_div(self, p, q, eps=1e-16):
+        return tnp.sum(p * tnp.log(p / q + eps))
+
     def isfinite(self, a):
         return tnp.isfinite(a)
 

ot/bregman.py

@@ -207,6 +207,9 @@ def sinkhorn2(a, b, M, reg, method='sinkhorn', numItermax=1000,
     - :math:`\mathbf{a}` and :math:`\mathbf{b}` are source and target
       weights (histograms, both sum to 1)
 
+    and returns :math:`\langle \gamma^*, \mathbf{M} \rangle_F` (without
+    the entropic contribution).
+
     .. note:: This function is backend-compatible and will work on arrays
         from all compatible backends.
 
@@ -320,15 +323,18 @@ def sinkhorn2(a, b, M, reg, method='sinkhorn', numItermax=1000,
     if len(b.shape) < 2:
         if method.lower() == 'sinkhorn':
             res = sinkhorn_knopp(a, b, M, reg, numItermax=numItermax,
-                                 stopThr=stopThr, verbose=verbose, log=log,
+                                 stopThr=stopThr, verbose=verbose,
+                                 log=log, warn=warn,
                                  **kwargs)
         elif method.lower() == 'sinkhorn_log':
             res = sinkhorn_log(a, b, M, reg, numItermax=numItermax,
-                               stopThr=stopThr, verbose=verbose, log=log,
+                               stopThr=stopThr, verbose=verbose,
+                               log=log, warn=warn,
                                **kwargs)
         elif method.lower() == 'sinkhorn_stabilized':
             res = sinkhorn_stabilized(a, b, M, reg, numItermax=numItermax,
-                                      stopThr=stopThr, verbose=verbose, log=log,
+                                      stopThr=stopThr, verbose=verbose,
+                                      log=log, warn=warn,
                                       **kwargs)
         else:
             raise ValueError("Unknown method '%s'." % method)
@@ -341,15 +347,18 @@ def sinkhorn2(a, b, M, reg, method='sinkhorn', numItermax=1000,
 
         if method.lower() == 'sinkhorn':
             return sinkhorn_knopp(a, b, M, reg, numItermax=numItermax,
-                                  stopThr=stopThr, verbose=verbose, log=log,
+                                  stopThr=stopThr, verbose=verbose,
+                                  log=log, warn=warn,
                                   **kwargs)
         elif method.lower() == 'sinkhorn_log':
             return sinkhorn_log(a, b, M, reg, numItermax=numItermax,
-                                stopThr=stopThr, verbose=verbose, log=log,
+                                stopThr=stopThr, verbose=verbose,
+                                log=log, warn=warn,
                                 **kwargs)
         elif method.lower() == 'sinkhorn_stabilized':
             return sinkhorn_stabilized(a, b, M, reg, numItermax=numItermax,
-                                       stopThr=stopThr, verbose=verbose, log=log,
+                                       stopThr=stopThr, verbose=verbose,
+                                       log=log, warn=warn,
                                        **kwargs)
         else:
             raise ValueError("Unknown method '%s'." % method)
@@ -1278,7 +1287,7 @@ def get_reg(n):  # exponential decreasing
         regi = get_reg(ii)
 
         G, logi = sinkhorn_stabilized(a, b, M, regi,
-                                      numItermax=numInnerItermax, stopThr=1e-9,
+                                      numItermax=numInnerItermax, stopThr=stopThr,
                                       warmstart=(alpha, beta), verbose=False,
                                       print_period=20, tau=tau, log=True)
 
@@ -3059,6 +3068,9 @@ def empirical_sinkhorn2(X_s, X_t, reg, a=None, b=None, metric='sqeuclidean',
       :math:`\Omega(\gamma)=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})`
     - :math:`\mathbf{a}` and :math:`\mathbf{b}` are source and target weights (sum to 1)
 
+    and returns :math:`\langle \gamma^*, \mathbf{M} \rangle_F` (without
+    the entropic contribution).
+
 
     Parameters
     ----------
@@ -3237,6 +3249,13 @@ def empirical_sinkhorn_divergence(X_s, X_t, reg, a=None, b=None, metric='sqeucli
       :math:`\Omega(\gamma)=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})`
     - :math:`\mathbf{a}` and :math:`\mathbf{b}` are source and target weights (sum to 1)
 
+    and returns :math:`\langle \gamma^*, \mathbf{M} \rangle_F -(\langle \gamma^*_a, \mathbf{M_a} \rangle_F + \langle
+    \gamma^*_b , \mathbf{M_b} \rangle_F)/2`.
+
+    .. note: The current implementation does not account for the entropic contributions and thus differs from the
+    Sinkhorn divergence as introduced in the literature. The possibility to account for the entropic contributions
+    will be provided in a future release.
+
 
     Parameters
     ----------
@@ -3293,17 +3312,17 @@ def empirical_sinkhorn_divergence(X_s, X_t, reg, a=None, b=None, metric='sqeucli
     if log:
         sinkhorn_loss_ab, log_ab = empirical_sinkhorn2(X_s, X_t, reg, a, b, metric=metric,
                                                        numIterMax=numIterMax,
-                                                       stopThr=1e-9, verbose=verbose,
+                                                       stopThr=stopThr, verbose=verbose,
                                                        log=log, warn=warn, **kwargs)
 
         sinkhorn_loss_a, log_a = empirical_sinkhorn2(X_s, X_s, reg, a, a, metric=metric,
                                                      numIterMax=numIterMax,
-                                                     stopThr=1e-9, verbose=verbose,
+                                                     stopThr=stopThr, verbose=verbose,
                                                      log=log, warn=warn, **kwargs)
 
         sinkhorn_loss_b, log_b = empirical_sinkhorn2(X_t, X_t, reg, b, b, metric=metric,
                                                      numIterMax=numIterMax,
-                                                     stopThr=1e-9, verbose=verbose,
+                                                     stopThr=stopThr, verbose=verbose,
                                                      log=log, warn=warn, **kwargs)
 
         sinkhorn_div = sinkhorn_loss_ab - 0.5 * (sinkhorn_loss_a + sinkhorn_loss_b)
@@ -3320,17 +3339,17 @@ def empirical_sinkhorn_divergence(X_s, X_t, reg, a=None, b=None, metric='sqeucli
 
     else:
         sinkhorn_loss_ab = empirical_sinkhorn2(X_s, X_t, reg, a, b, metric=metric,
-                                               numIterMax=numIterMax, stopThr=1e-9,
+                                               numIterMax=numIterMax, stopThr=stopThr,
                                                verbose=verbose, log=log,
                                                warn=warn, **kwargs)
 
         sinkhorn_loss_a = empirical_sinkhorn2(X_s, X_s, reg, a, a, metric=metric,
-                                              numIterMax=numIterMax, stopThr=1e-9,
+                                              numIterMax=numIterMax, stopThr=stopThr,
                                               verbose=verbose, log=log,
                                               warn=warn, **kwargs)
 
         sinkhorn_loss_b = empirical_sinkhorn2(X_t, X_t, reg, b, b, metric=metric,
-                                              numIterMax=numIterMax, stopThr=1e-9,
+                                              numIterMax=numIterMax, stopThr=stopThr,
                                               verbose=verbose, log=log,
                                               warn=warn, **kwargs)
 

ot/da.py

@@ -126,8 +126,12 @@ def sinkhorn_lpl1_mm(a, labels_a, b, M, reg, eta=0.1, numItermax=10,
     W = nx.zeros(M.shape, type_as=M)
     for cpt in range(numItermax):
         Mreg = M + eta * W
-        transp = sinkhorn(a, b, Mreg, reg, numItermax=numInnerItermax,
-                          stopThr=stopInnerThr)
+        if log:
+            transp, log = sinkhorn(a, b, Mreg, reg, numItermax=numInnerItermax,
+                                   stopThr=stopInnerThr, log=True)
+        else:
+            transp = sinkhorn(a, b, Mreg, reg, numItermax=numInnerItermax,
+                              stopThr=stopInnerThr)
         # the transport has been computed. Check if classes are really
         # separated
         W = nx.ones(M.shape, type_as=M)
@@ -136,7 +140,10 @@ def sinkhorn_lpl1_mm(a, labels_a, b, M, reg, eta=0.1, numItermax=10,
             majs = p * ((majs + epsilon) ** (p - 1))
             W[indices_labels[i]] = majs
 
-    return transp
+    if log:
+        return transp, log
+    else:
+        return transp
 
 
 def sinkhorn_l1l2_gl(a, labels_a, b, M, reg, eta=0.1, numItermax=10,

ot/helpers/pre_build_helpers.py

@@ -4,34 +4,14 @@
 import sys
 import glob
 import tempfile
-import setuptools  # noqa
 import subprocess
 
-from distutils.dist import Distribution
-from distutils.sysconfig import customize_compiler
-from numpy.distutils.ccompiler import new_compiler
-from numpy.distutils.command.config_compiler import config_cc
+from setuptools.command.build_ext import customize_compiler, new_compiler
 
 
 def _get_compiler():
-    """Get a compiler equivalent to the one that will be used to build POT
-    Handles compiler specified as follows:
-        - python setup.py build_ext --compiler=<compiler>
-        - CC=<compiler> python setup.py build_ext
-    """
-    dist = Distribution({'script_name': os.path.basename(sys.argv[0]),
-                         'script_args': sys.argv[1:],
-                         'cmdclass': {'config_cc': config_cc}})
-
-    cmd_opts = dist.command_options.get('build_ext')
-    if cmd_opts is not None and 'compiler' in cmd_opts:
-        compiler = cmd_opts['compiler'][1]
-    else:
-        compiler = None
-
-    ccompiler = new_compiler(compiler=compiler)
+    ccompiler = new_compiler()
     customize_compiler(ccompiler)
-
     return ccompiler
 
 

ot/lp/network_simplex_simple_omp.h

@@ -67,7 +67,7 @@
 //#include "core.h"
 //#include "lmath.h"
 
-#ifdef OMP
+#ifdef _OPENMP
 #include <omp.h>
 #endif
 #include <cmath>
@@ -254,7 +254,7 @@ namespace lemon_omp {
 			// Reset data structures
 			reset();
 			max_iter = maxiters;
-#ifdef OMP
+#ifdef _OPENMP
             if (max_threads < 0) {
                 max_threads = omp_get_max_threads();
             }
@@ -513,7 +513,7 @@ namespace lemon_omp {
 					int j;
 #pragma omp parallel
 					{
-#ifdef OMP
+#ifdef _OPENMP
 						int t = omp_get_thread_num();
 #else
 						int t = 0;