Merge branch 'master' into doc_modules

6 files changed, 191 insertions, 172 deletions

diff --git a/ot/bregman.py b/ot/bregman.py
index 09716e6..13dfa3b 100644
--- a/ot/bregman.py
+++ b/ot/bregman.py
@@ -16,7 +16,7 @@ from .utils import unif, dist
 
 def sinkhorn(a, b, M, reg, method='sinkhorn', numItermax=1000,
              stopThr=1e-9, verbose=False, log=False, **kwargs):
-    u"""
+    r"""
     Solve the entropic regularization optimal transport problem and return the OT matrix
 
     The function solves the following optimization problem:
@@ -73,12 +73,12 @@ def sinkhorn(a, b, M, reg, method='sinkhorn', numItermax=1000,
     --------
 
     >>> import ot
-    >>> a=[.5,.5]
-    >>> b=[.5,.5]
-    >>> M=[[0.,1.],[1.,0.]]
-    >>> ot.sinkhorn(a,b,M,1)
-    array([[ 0.36552929,  0.13447071],
-           [ 0.13447071,  0.36552929]])
+    >>> a=[.5, .5]
+    >>> b=[.5, .5]
+    >>> M=[[0., 1.], [1., 0.]]
+    >>> ot.sinkhorn(a, b, M, 1)
+    array([[0.36552929, 0.13447071],
+           [0.13447071, 0.36552929]])
 
 
     References
@@ -131,7 +131,7 @@ def sinkhorn(a, b, M, reg, method='sinkhorn', numItermax=1000,
 
 def sinkhorn2(a, b, M, reg, method='sinkhorn', numItermax=1000,
               stopThr=1e-9, verbose=False, log=False, **kwargs):
-    u"""
+    r"""
     Solve the entropic regularization optimal transport problem and return the loss
 
     The function solves the following optimization problem:
@@ -188,11 +188,11 @@ def sinkhorn2(a, b, M, reg, method='sinkhorn', numItermax=1000,
     --------
 
     >>> import ot
-    >>> a=[.5,.5]
-    >>> b=[.5,.5]
-    >>> M=[[0.,1.],[1.,0.]]
-    >>> ot.sinkhorn2(a,b,M,1)
-    array([ 0.26894142])
+    >>> a=[.5, .5]
+    >>> b=[.5, .5]
+    >>> M=[[0., 1.], [1., 0.]]
+    >>> ot.sinkhorn2(a, b, M, 1)
+    array([0.26894142])
 
 
 
@@ -248,7 +248,7 @@ def sinkhorn2(a, b, M, reg, method='sinkhorn', numItermax=1000,
 
 def sinkhorn_knopp(a, b, M, reg, numItermax=1000,
                    stopThr=1e-9, verbose=False, log=False, **kwargs):
-    """
+    r"""
     Solve the entropic regularization optimal transport problem and return the OT matrix
 
     The function solves the following optimization problem:
@@ -302,12 +302,12 @@ def sinkhorn_knopp(a, b, M, reg, numItermax=1000,
     --------
 
     >>> import ot
-    >>> a=[.5,.5]
-    >>> b=[.5,.5]
-    >>> M=[[0.,1.],[1.,0.]]
-    >>> ot.sinkhorn(a,b,M,1)
-    array([[ 0.36552929,  0.13447071],
-           [ 0.13447071,  0.36552929]])
+    >>> a=[.5, .5]
+    >>> b=[.5, .5]
+    >>> M=[[0., 1.], [1., 0.]]
+    >>> ot.sinkhorn(a, b, M, 1)
+    array([[0.36552929, 0.13447071],
+           [0.13447071, 0.36552929]])
 
 
     References
@@ -422,7 +422,7 @@ def sinkhorn_knopp(a, b, M, reg, numItermax=1000,
 
 
 def greenkhorn(a, b, M, reg, numItermax=10000, stopThr=1e-9, verbose=False, log=False):
-    """
+    r"""
     Solve the entropic regularization optimal transport problem and return the OT matrix
 
     The algorithm used is based on the paper
@@ -481,12 +481,12 @@ def greenkhorn(a, b, M, reg, numItermax=10000, stopThr=1e-9, verbose=False, log=
     --------
 
     >>> import ot
-    >>> a=[.5,.5]
-    >>> b=[.5,.5]
-    >>> M=[[0.,1.],[1.,0.]]
-    >>> ot.bregman.greenkhorn(a,b,M,1)
-    array([[ 0.36552929,  0.13447071],
-           [ 0.13447071,  0.36552929]])
+    >>> a=[.5, .5]
+    >>> b=[.5, .5]
+    >>> M=[[0., 1.], [1., 0.]]
+    >>> ot.bregman.greenkhorn(a, b, M, 1)
+    array([[0.36552929, 0.13447071],
+           [0.13447071, 0.36552929]])
 
 
     References
@@ -576,7 +576,7 @@ def greenkhorn(a, b, M, reg, numItermax=10000, stopThr=1e-9, verbose=False, log=
 
 def sinkhorn_stabilized(a, b, M, reg, numItermax=1000, tau=1e3, stopThr=1e-9,
                         warmstart=None, verbose=False, print_period=20, log=False, **kwargs):
-    """
+    r"""
     Solve the entropic regularization OT problem with log stabilization
 
     The function solves the following optimization problem:
@@ -639,8 +639,8 @@ def sinkhorn_stabilized(a, b, M, reg, numItermax=1000, tau=1e3, stopThr=1e-9,
     >>> b=[.5,.5]
     >>> M=[[0.,1.],[1.,0.]]
     >>> ot.bregman.sinkhorn_stabilized(a,b,M,1)
-    array([[ 0.36552929,  0.13447071],
-           [ 0.13447071,  0.36552929]])
+    array([[0.36552929, 0.13447071],
+           [0.13447071, 0.36552929]])
 
 
     References
@@ -796,7 +796,7 @@ def sinkhorn_stabilized(a, b, M, reg, numItermax=1000, tau=1e3, stopThr=1e-9,
 
 def sinkhorn_epsilon_scaling(a, b, M, reg, numItermax=100, epsilon0=1e4, numInnerItermax=100,
                              tau=1e3, stopThr=1e-9, warmstart=None, verbose=False, print_period=10, log=False, **kwargs):
-    """
+    r"""
     Solve the entropic regularization optimal transport problem with log
     stabilization and epsilon scaling.
 
@@ -862,12 +862,12 @@ def sinkhorn_epsilon_scaling(a, b, M, reg, numItermax=100, epsilon0=1e4, numInne
     --------
 
     >>> import ot
-    >>> a=[.5,.5]
-    >>> b=[.5,.5]
-    >>> M=[[0.,1.],[1.,0.]]
-    >>> ot.bregman.sinkhorn_epsilon_scaling(a,b,M,1)
-    array([[ 0.36552929,  0.13447071],
-           [ 0.13447071,  0.36552929]])
+    >>> a=[.5, .5]
+    >>> b=[.5, .5]
+    >>> M=[[0., 1.], [1., 0.]]
+    >>> ot.bregman.sinkhorn_epsilon_scaling(a, b, M, 1)
+    array([[0.36552929, 0.13447071],
+           [0.13447071, 0.36552929]])
 
 
     References
@@ -989,7 +989,7 @@ def projC(gamma, q):
 
 def barycenter(A, M, reg, weights=None, numItermax=1000,
                stopThr=1e-4, verbose=False, log=False):
-    """Compute the entropic regularized wasserstein barycenter of distributions A
+    r"""Compute the entropic regularized wasserstein barycenter of distributions A
 
      The function solves the following optimization problem:
 
@@ -1084,7 +1084,7 @@ def barycenter(A, M, reg, weights=None, numItermax=1000,
 
 
 def convolutional_barycenter2d(A, reg, weights=None, numItermax=10000, stopThr=1e-9, stabThr=1e-30, verbose=False, log=False):
-    """Compute the entropic regularized wasserstein barycenter of distributions A
+    r"""Compute the entropic regularized wasserstein barycenter of distributions A
     where A is a collection of 2D images.
 
      The function solves the following optimization problem:
@@ -1195,7 +1195,7 @@ def convolutional_barycenter2d(A, reg, weights=None, numItermax=10000, stopThr=1
 
 def unmix(a, D, M, M0, h0, reg, reg0, alpha, numItermax=1000,
           stopThr=1e-3, verbose=False, log=False):
-    """
+    r"""
     Compute the unmixing of an observation with a given dictionary using Wasserstein distance
 
     The function solve the following optimization problem:
@@ -1302,7 +1302,7 @@ def unmix(a, D, M, M0, h0, reg, reg0, alpha, numItermax=1000,
 
 
 def empirical_sinkhorn(X_s, X_t, reg, a=None, b=None, metric='sqeuclidean', numIterMax=10000, stopThr=1e-9, verbose=False, log=False, **kwargs):
-    '''
+    r'''
     Solve the entropic regularization optimal transport problem and return the
     OT matrix from empirical data
 
@@ -1360,10 +1360,9 @@ def empirical_sinkhorn(X_s, X_t, reg, a=None, b=None, metric='sqeuclidean', numI
     >>> reg = 0.1
     >>> X_s = np.reshape(np.arange(n_s), (n_s, 1))
     >>> X_t = np.reshape(np.arange(0, n_t), (n_t, 1))
-    >>> emp_sinkhorn = empirical_sinkhorn(X_s, X_t, reg, verbose=False)
-    >>> print(emp_sinkhorn)
-    >>> [[4.99977301e-01 2.26989344e-05]
-        [2.26989344e-05 4.99977301e-01]]
+    >>> empirical_sinkhorn(X_s, X_t, reg, verbose=False)  # doctest: +NORMALIZE_WHITESPACE
+    array([[4.99977301e-01,  2.26989344e-05],
+           [2.26989344e-05,  4.99977301e-01]])
 
 
     References
@@ -1392,7 +1391,7 @@ def empirical_sinkhorn(X_s, X_t, reg, a=None, b=None, metric='sqeuclidean', numI
 
 
 def empirical_sinkhorn2(X_s, X_t, reg, a=None, b=None, metric='sqeuclidean', numIterMax=10000, stopThr=1e-9, verbose=False, log=False, **kwargs):
-    '''
+    r'''
     Solve the entropic regularization optimal transport problem from empirical
     data and return the OT loss
 
@@ -1451,9 +1450,8 @@ def empirical_sinkhorn2(X_s, X_t, reg, a=None, b=None, metric='sqeuclidean', num
     >>> reg = 0.1
     >>> X_s = np.reshape(np.arange(n_s), (n_s, 1))
     >>> X_t = np.reshape(np.arange(0, n_t), (n_t, 1))
-    >>> loss_sinkhorn = empirical_sinkhorn2(X_s, X_t, reg, verbose=False)
-    >>> print(loss_sinkhorn)
-    >>> [4.53978687e-05]
+    >>> empirical_sinkhorn2(X_s, X_t, reg, verbose=False)
+    array([4.53978687e-05])
 
 
     References
@@ -1482,7 +1480,7 @@ def empirical_sinkhorn2(X_s, X_t, reg, a=None, b=None, metric='sqeuclidean', num
 
 
 def empirical_sinkhorn_divergence(X_s, X_t, reg, a=None, b=None, metric='sqeuclidean', numIterMax=10000, stopThr=1e-9, verbose=False, log=False, **kwargs):
-    '''
+    r'''
     Compute the sinkhorn divergence loss from empirical data
 
     The function solves the following optimization problems and return the
@@ -1560,9 +1558,8 @@ def empirical_sinkhorn_divergence(X_s, X_t, reg, a=None, b=None, metric='sqeucli
     >>> reg = 0.1
     >>> X_s = np.reshape(np.arange(n_s), (n_s, 1))
     >>> X_t = np.reshape(np.arange(0, n_t), (n_t, 1))
-    >>> emp_sinkhorn_div = empirical_sinkhorn_divergence(X_s, X_t, reg)
-    >>> print(emp_sinkhorn_div)
-    >>> [2.99977435]
+    >>> empirical_sinkhorn_divergence(X_s, X_t, reg)  # doctest: +ELLIPSIS
+    array([1.499...])
 
 
     References
diff --git a/ot/gpu/bregman.py b/ot/gpu/bregman.py
index 978b307..2e2df83 100644
--- a/ot/gpu/bregman.py
+++ b/ot/gpu/bregman.py
@@ -70,17 +70,6 @@ def sinkhorn_knopp(a, b, M, reg, numItermax=1000, stopThr=1e-9,
     log : dict
         log dictionary return only if log==True in parameters
 
-    Examples
-    --------
-
-    >>> import ot
-    >>> a=[.5,.5]
-    >>> b=[.5,.5]
-    >>> M=[[0.,1.],[1.,0.]]
-    >>> ot.sinkhorn(a,b,M,1)
-    array([[ 0.36552929,  0.13447071],
-           [ 0.13447071,  0.36552929]])
-
 
     References
     ----------
diff --git a/ot/lp/__init__.py b/ot/lp/__init__.py
index 8bbd6d9..aed29f8 100644
--- a/ot/lp/__init__.py
+++ b/ot/lp/__init__.py
@@ -28,7 +28,7 @@ __all__=['emd', 'emd2', 'barycenter', 'free_support_barycenter', 'cvx',
 
 
 def emd(a, b, M, numItermax=100000, log=False):
-    """Solves the Earth Movers distance problem and returns the OT matrix
+    r"""Solves the Earth Movers distance problem and returns the OT matrix
 
 
     .. math::
@@ -83,8 +83,8 @@ def emd(a, b, M, numItermax=100000, log=False):
     >>> b=[.5,.5]
     >>> M=[[0.,1.],[1.,0.]]
     >>> ot.emd(a,b,M)
-    array([[ 0.5,  0. ],
-           [ 0. ,  0.5]])
+    array([[0.5, 0. ],
+           [0. , 0.5]])
 
     References
     ----------
@@ -124,7 +124,7 @@ def emd(a, b, M, numItermax=100000, log=False):
 
 def emd2(a, b, M, processes=multiprocessing.cpu_count(),
          numItermax=100000, log=False, return_matrix=False):
-    """Solves the Earth Movers distance problem and returns the loss
+    r"""Solves the Earth Movers distance problem and returns the loss
 
     .. math::
         \gamma = arg\min_\gamma <\gamma,M>_F
@@ -326,7 +326,7 @@ def free_support_barycenter(measures_locations, measures_weights, X_init, b=None
 
 def emd_1d(x_a, x_b, a=None, b=None, metric='sqeuclidean', p=1., dense=True,
            log=False):
-    """Solves the Earth Movers distance problem between 1d measures and returns
+    r"""Solves the Earth Movers distance problem between 1d measures and returns
     the OT matrix
 
 
@@ -392,11 +392,11 @@ def emd_1d(x_a, x_b, a=None, b=None, metric='sqeuclidean', p=1., dense=True,
     >>> x_a = [2., 0.]
     >>> x_b = [0., 3.]
     >>> ot.emd_1d(x_a, x_b, a, b)
-    array([[0. ,  0.5],
-           [0.5,  0. ]])
+    array([[0. , 0.5],
+           [0.5, 0. ]])
     >>> ot.emd_1d(x_a, x_b)
-    array([[0. ,  0.5],
-           [0.5,  0. ]])
+    array([[0. , 0.5],
+           [0.5, 0. ]])
 
     References
     ----------
@@ -446,7 +446,7 @@ def emd_1d(x_a, x_b, a=None, b=None, metric='sqeuclidean', p=1., dense=True,
 
 def emd2_1d(x_a, x_b, a=None, b=None, metric='sqeuclidean', p=1., dense=True,
             log=False):
-    """Solves the Earth Movers distance problem between 1d measures and returns
+    r"""Solves the Earth Movers distance problem between 1d measures and returns
     the loss
 
 
@@ -541,7 +541,7 @@ 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.):
-    """Solves the p-Wasserstein distance problem between 1d measures and returns
+    r"""Solves the p-Wasserstein distance problem between 1d measures and returns
     the distance
 
 
diff --git a/ot/stochastic.py b/ot/stochastic.py
index 85c4230..bf3e7a7 100644
--- a/ot/stochastic.py
+++ b/ot/stochastic.py
@@ -11,7 +11,7 @@ import numpy as np
 
 
 def coordinate_grad_semi_dual(b, M, reg, beta, i):
-    '''
+    r'''
     Compute the coordinate gradient update for regularized discrete distributions for (i, :)
 
     The function computes the gradient of the semi dual problem:
@@ -51,21 +51,24 @@ def coordinate_grad_semi_dual(b, M, reg, beta, i):
 
     Examples
     --------
-
+    >>> import ot
+    >>> np.random.seed(0)
     >>> n_source = 7
     >>> n_target = 4
-    >>> reg = 1
-    >>> numItermax = 300000
     >>> a = ot.utils.unif(n_source)
     >>> b = ot.utils.unif(n_target)
-    >>> rng = np.random.RandomState(0)
-    >>> X_source = rng.randn(n_source, 2)
-    >>> Y_target = rng.randn(n_target, 2)
+    >>> X_source = np.random.randn(n_source, 2)
+    >>> Y_target = np.random.randn(n_target, 2)
     >>> M = ot.dist(X_source, Y_target)
-    >>> method = "ASGD"
-    >>> asgd_pi = stochastic.solve_semi_dual_entropic(a, b, M, reg,
-                                                      method, numItermax)
-    >>> print(asgd_pi)
+    >>> ot.stochastic.solve_semi_dual_entropic(a, b, M, reg=1, method="ASGD", numItermax=300000)
+    array([[2.53942342e-02, 9.98640673e-02, 1.75945647e-02, 4.27664307e-06],
+           [1.21556999e-01, 1.26350515e-02, 1.30491795e-03, 7.36017394e-03],
+           [3.54070702e-03, 7.63581358e-02, 6.29581672e-02, 1.32812798e-07],
+           [2.60578198e-02, 3.35916645e-02, 8.28023223e-02, 4.05336238e-04],
+           [9.86808864e-03, 7.59774324e-04, 1.08702729e-02, 1.21359007e-01],
+           [2.17218856e-02, 9.12931802e-04, 1.87962526e-03, 1.18342700e-01],
+           [4.14237512e-02, 2.67487857e-02, 7.23016955e-02, 2.38291052e-03]])
+
 
     References
     ----------
@@ -84,7 +87,7 @@ def coordinate_grad_semi_dual(b, M, reg, beta, i):
 
 
 def sag_entropic_transport(a, b, M, reg, numItermax=10000, lr=None):
-    '''
+    r'''
     Compute the SAG algorithm to solve the regularized discrete measures
         optimal transport max problem
 
@@ -133,21 +136,23 @@ def sag_entropic_transport(a, b, M, reg, numItermax=10000, lr=None):
 
     Examples
     --------
-
+    >>> import ot
+    >>> np.random.seed(0)
     >>> n_source = 7
     >>> n_target = 4
-    >>> reg = 1
-    >>> numItermax = 300000
     >>> a = ot.utils.unif(n_source)
     >>> b = ot.utils.unif(n_target)
-    >>> rng = np.random.RandomState(0)
-    >>> X_source = rng.randn(n_source, 2)
-    >>> Y_target = rng.randn(n_target, 2)
+    >>> X_source = np.random.randn(n_source, 2)
+    >>> Y_target = np.random.randn(n_target, 2)
     >>> M = ot.dist(X_source, Y_target)
-    >>> method = "ASGD"
-    >>> asgd_pi = stochastic.solve_semi_dual_entropic(a, b, M, reg,
-                                                      method, numItermax)
-    >>> print(asgd_pi)
+    >>> ot.stochastic.solve_semi_dual_entropic(a, b, M, reg=1, method="ASGD", numItermax=300000)
+    array([[2.53942342e-02, 9.98640673e-02, 1.75945647e-02, 4.27664307e-06],
+           [1.21556999e-01, 1.26350515e-02, 1.30491795e-03, 7.36017394e-03],
+           [3.54070702e-03, 7.63581358e-02, 6.29581672e-02, 1.32812798e-07],
+           [2.60578198e-02, 3.35916645e-02, 8.28023223e-02, 4.05336238e-04],
+           [9.86808864e-03, 7.59774324e-04, 1.08702729e-02, 1.21359007e-01],
+           [2.17218856e-02, 9.12931802e-04, 1.87962526e-03, 1.18342700e-01],
+           [4.14237512e-02, 2.67487857e-02, 7.23016955e-02, 2.38291052e-03]])
 
     References
     ----------
@@ -176,7 +181,7 @@ def sag_entropic_transport(a, b, M, reg, numItermax=10000, lr=None):
 
 
 def averaged_sgd_entropic_transport(a, b, M, reg, numItermax=300000, lr=None):
-    '''
+    r'''
     Compute the ASGD algorithm to solve the regularized semi continous measures optimal transport max problem
 
     The function solves the following optimization problem:
@@ -223,21 +228,23 @@ def averaged_sgd_entropic_transport(a, b, M, reg, numItermax=300000, lr=None):
 
     Examples
     --------
-
+    >>> import ot
+    >>> np.random.seed(0)
     >>> n_source = 7
     >>> n_target = 4
-    >>> reg = 1
-    >>> numItermax = 300000
     >>> a = ot.utils.unif(n_source)
     >>> b = ot.utils.unif(n_target)
-    >>> rng = np.random.RandomState(0)
-    >>> X_source = rng.randn(n_source, 2)
-    >>> Y_target = rng.randn(n_target, 2)
+    >>> X_source = np.random.randn(n_source, 2)
+    >>> Y_target = np.random.randn(n_target, 2)
     >>> M = ot.dist(X_source, Y_target)
-    >>> method = "ASGD"
-    >>> asgd_pi = stochastic.solve_semi_dual_entropic(a, b, M, reg,
-                                                      method, numItermax)
-    >>> print(asgd_pi)
+    >>> ot.stochastic.solve_semi_dual_entropic(a, b, M, reg=1, method="ASGD", numItermax=300000)
+    array([[2.53942342e-02, 9.98640673e-02, 1.75945647e-02, 4.27664307e-06],
+           [1.21556999e-01, 1.26350515e-02, 1.30491795e-03, 7.36017394e-03],
+           [3.54070702e-03, 7.63581358e-02, 6.29581672e-02, 1.32812798e-07],
+           [2.60578198e-02, 3.35916645e-02, 8.28023223e-02, 4.05336238e-04],
+           [9.86808864e-03, 7.59774324e-04, 1.08702729e-02, 1.21359007e-01],
+           [2.17218856e-02, 9.12931802e-04, 1.87962526e-03, 1.18342700e-01],
+           [4.14237512e-02, 2.67487857e-02, 7.23016955e-02, 2.38291052e-03]])
 
     References
     ----------
@@ -264,7 +271,7 @@ def averaged_sgd_entropic_transport(a, b, M, reg, numItermax=300000, lr=None):
 
 
 def c_transform_entropic(b, M, reg, beta):
-    '''
+    r'''
     The goal is to recover u from the c-transform.
 
     The function computes the c_transform of a dual variable from the other
@@ -303,21 +310,23 @@ def c_transform_entropic(b, M, reg, beta):
 
     Examples
     --------
-
+    >>> import ot
+    >>> np.random.seed(0)
     >>> n_source = 7
     >>> n_target = 4
-    >>> reg = 1
-    >>> numItermax = 300000
     >>> a = ot.utils.unif(n_source)
     >>> b = ot.utils.unif(n_target)
-    >>> rng = np.random.RandomState(0)
-    >>> X_source = rng.randn(n_source, 2)
-    >>> Y_target = rng.randn(n_target, 2)
+    >>> X_source = np.random.randn(n_source, 2)
+    >>> Y_target = np.random.randn(n_target, 2)
     >>> M = ot.dist(X_source, Y_target)
-    >>> method = "ASGD"
-    >>> asgd_pi = stochastic.solve_semi_dual_entropic(a, b, M, reg,
-                                                      method, numItermax)
-    >>> print(asgd_pi)
+    >>> ot.stochastic.solve_semi_dual_entropic(a, b, M, reg=1, method="ASGD", numItermax=300000)
+    array([[2.53942342e-02, 9.98640673e-02, 1.75945647e-02, 4.27664307e-06],
+           [1.21556999e-01, 1.26350515e-02, 1.30491795e-03, 7.36017394e-03],
+           [3.54070702e-03, 7.63581358e-02, 6.29581672e-02, 1.32812798e-07],
+           [2.60578198e-02, 3.35916645e-02, 8.28023223e-02, 4.05336238e-04],
+           [9.86808864e-03, 7.59774324e-04, 1.08702729e-02, 1.21359007e-01],
+           [2.17218856e-02, 9.12931802e-04, 1.87962526e-03, 1.18342700e-01],
+           [4.14237512e-02, 2.67487857e-02, 7.23016955e-02, 2.38291052e-03]])
 
     References
     ----------
@@ -340,7 +349,7 @@ def c_transform_entropic(b, M, reg, beta):
 
 def solve_semi_dual_entropic(a, b, M, reg, method, numItermax=10000, lr=None,
                                 log=False):
-    '''
+    r'''
     Compute the transportation matrix to solve the regularized discrete
         measures optimal transport max problem
 
@@ -398,21 +407,23 @@ def solve_semi_dual_entropic(a, b, M, reg, method, numItermax=10000, lr=None,
 
     Examples
     --------
-
+    >>> import ot
+    >>> np.random.seed(0)
     >>> n_source = 7
     >>> n_target = 4
-    >>> reg = 1
-    >>> numItermax = 300000
     >>> a = ot.utils.unif(n_source)
     >>> b = ot.utils.unif(n_target)
-    >>> rng = np.random.RandomState(0)
-    >>> X_source = rng.randn(n_source, 2)
-    >>> Y_target = rng.randn(n_target, 2)
+    >>> X_source = np.random.randn(n_source, 2)
+    >>> Y_target = np.random.randn(n_target, 2)
     >>> M = ot.dist(X_source, Y_target)
-    >>> method = "ASGD"
-    >>> asgd_pi = stochastic.solve_semi_dual_entropic(a, b, M, reg,
-                                                      method, numItermax)
-    >>> print(asgd_pi)
+    >>> ot.stochastic.solve_semi_dual_entropic(a, b, M, reg=1, method="ASGD", numItermax=300000)
+    array([[2.53942342e-02, 9.98640673e-02, 1.75945647e-02, 4.27664307e-06],
+           [1.21556999e-01, 1.26350515e-02, 1.30491795e-03, 7.36017394e-03],
+           [3.54070702e-03, 7.63581358e-02, 6.29581672e-02, 1.32812798e-07],
+           [2.60578198e-02, 3.35916645e-02, 8.28023223e-02, 4.05336238e-04],
+           [9.86808864e-03, 7.59774324e-04, 1.08702729e-02, 1.21359007e-01],
+           [2.17218856e-02, 9.12931802e-04, 1.87962526e-03, 1.18342700e-01],
+           [4.14237512e-02, 2.67487857e-02, 7.23016955e-02, 2.38291052e-03]])
 
     References
     ----------
@@ -451,7 +462,7 @@ def solve_semi_dual_entropic(a, b, M, reg, method, numItermax=10000, lr=None,
 
 def batch_grad_dual(a, b, M, reg, alpha, beta, batch_size, batch_alpha,
                     batch_beta):
-    '''
+    r'''
     Computes the partial gradient of the dual optimal transport problem.
 
     For each (i,j) in a batch of coordinates, the partial gradients are :
@@ -506,25 +517,29 @@ def batch_grad_dual(a, b, M, reg, alpha, beta, batch_size, batch_alpha,
 
     Examples
     --------
-
+    >>> import ot
+    >>> np.random.seed(0)
     >>> n_source = 7
     >>> n_target = 4
-    >>> reg = 1
-    >>> numItermax = 20000
-    >>> lr = 0.1
-    >>> batch_size = 3
-    >>> log = True
     >>> a = ot.utils.unif(n_source)
     >>> b = ot.utils.unif(n_target)
-    >>> rng = np.random.RandomState(0)
-    >>> X_source = rng.randn(n_source, 2)
-    >>> Y_target = rng.randn(n_target, 2)
+    >>> X_source = np.random.randn(n_source, 2)
+    >>> Y_target = np.random.randn(n_target, 2)
     >>> M = ot.dist(X_source, Y_target)
-    >>> sgd_dual_pi, log = stochastic.solve_dual_entropic(a, b, M, reg,
-                                                            batch_size,
-                                                            numItermax, lr, log)
-    >>> print(log['alpha'], log['beta'])
-    >>> print(sgd_dual_pi)
+    >>> sgd_dual_pi, log = ot.stochastic.solve_dual_entropic(a, b, M, reg=1, batch_size=3, numItermax=30000, lr=0.1, log=True)
+    >>> log['alpha']
+    array([0.71759102, 1.57057384, 0.85576566, 0.1208211 , 0.59190466,
+           1.197148  , 0.17805133])
+    >>> log['beta']
+    array([0.49741367, 0.57478564, 1.40075528, 2.75890102])
+    >>> sgd_dual_pi
+    array([[2.09730063e-02, 8.38169324e-02, 7.50365455e-03, 8.72731415e-09],
+           [5.58432437e-03, 5.89881299e-04, 3.09558411e-05, 8.35469849e-07],
+           [3.26489515e-03, 7.15536035e-02, 2.99778211e-02, 3.02601593e-10],
+           [4.05390622e-02, 5.31085068e-02, 6.65191787e-02, 1.55812785e-06],
+           [7.82299812e-02, 6.12099102e-03, 4.44989098e-02, 2.37719187e-03],
+           [5.06266486e-02, 2.16230494e-03, 2.26215141e-03, 6.81514609e-04],
+           [6.06713990e-02, 3.98139808e-02, 5.46829338e-02, 8.62371424e-06]])
 
     References
     ----------
@@ -533,7 +548,6 @@ def batch_grad_dual(a, b, M, reg, alpha, beta, batch_size, batch_alpha,
                     International Conference on Learning Representation (2018),
                       arXiv preprint arxiv:1711.02283.
     '''
-
     G = - (np.exp((alpha[batch_alpha, None] + beta[None, batch_beta] -
                    M[batch_alpha, :][:, batch_beta]) / reg) *
            a[batch_alpha, None] * b[None, batch_beta])
@@ -548,7 +562,7 @@ def batch_grad_dual(a, b, M, reg, alpha, beta, batch_size, batch_alpha,
 
 
 def sgd_entropic_regularization(a, b, M, reg, batch_size, numItermax, lr):
-    '''
+    r'''
     Compute the sgd algorithm to solve the regularized discrete measures
         optimal transport dual problem
 
@@ -597,7 +611,7 @@ def sgd_entropic_regularization(a, b, M, reg, batch_size, numItermax, lr):
 
     Examples
     --------
-
+    >>> import ot
     >>> n_source = 7
     >>> n_target = 4
     >>> reg = 1
@@ -611,11 +625,20 @@ def sgd_entropic_regularization(a, b, M, reg, batch_size, numItermax, lr):
     >>> X_source = rng.randn(n_source, 2)
     >>> Y_target = rng.randn(n_target, 2)
     >>> M = ot.dist(X_source, Y_target)
-    >>> sgd_dual_pi, log = stochastic.solve_dual_entropic(a, b, M, reg,
-                                                          batch_size,
-                                                          numItermax, lr, log)
-    >>> print(log['alpha'], log['beta'])
-    >>> print(sgd_dual_pi)
+    >>> sgd_dual_pi, log = ot.stochastic.solve_dual_entropic(a, b, M, reg, batch_size, numItermax, lr, log)
+    >>> log['alpha']
+    array([0.64171798, 1.27932201, 0.78132257, 0.15638935, 0.54888354,
+           1.03663469, 0.20595781])
+    >>> log['beta']
+    array([0.51207194, 0.58033189, 1.28922676, 2.26859736])
+    >>> sgd_dual_pi
+    array([[1.97276541e-02, 7.81248547e-02, 6.22136048e-03, 4.95442423e-09],
+           [4.23494310e-03, 4.43286263e-04, 2.06927079e-05, 3.82389139e-07],
+           [3.07542414e-03, 6.67897769e-02, 2.48904999e-02, 1.72030247e-10],
+           [4.26271990e-02, 5.53375455e-02, 6.16535024e-02, 9.88812650e-07],
+           [7.60423265e-02, 5.89585256e-03, 3.81267087e-02, 1.39458256e-03],
+           [4.37557504e-02, 1.85189176e-03, 1.72335760e-03, 3.55491279e-04],
+           [6.33096109e-02, 4.11683954e-02, 5.02962051e-02, 5.43097516e-06]])
 
     References
     ----------
@@ -644,7 +667,7 @@ def sgd_entropic_regularization(a, b, M, reg, batch_size, numItermax, lr):
 
 def solve_dual_entropic(a, b, M, reg, batch_size, numItermax=10000, lr=1,
                         log=False):
-    '''
+    r'''
     Compute the transportation matrix to solve the regularized discrete measures
         optimal transport dual problem
 
@@ -695,7 +718,7 @@ def solve_dual_entropic(a, b, M, reg, batch_size, numItermax=10000, lr=1,
 
     Examples
     --------
-
+    >>> import ot
     >>> n_source = 7
     >>> n_target = 4
     >>> reg = 1
@@ -709,11 +732,20 @@ def solve_dual_entropic(a, b, M, reg, batch_size, numItermax=10000, lr=1,
     >>> X_source = rng.randn(n_source, 2)
     >>> Y_target = rng.randn(n_target, 2)
     >>> M = ot.dist(X_source, Y_target)
-    >>> sgd_dual_pi, log = stochastic.solve_dual_entropic(a, b, M, reg,
-                                                            batch_size,
-                                                            numItermax, lr, log)
-    >>> print(log['alpha'], log['beta'])
-    >>> print(sgd_dual_pi)
+    >>> sgd_dual_pi, log = ot.stochastic.solve_dual_entropic(a, b, M, reg, batch_size, numItermax, lr, log)
+    >>> log['alpha']
+    array([0.64057733, 1.2683513 , 0.75610161, 0.16024284, 0.54926534,
+           1.0514201 , 0.19958936])
+    >>> log['beta']
+    array([0.51372571, 0.58843489, 1.27993921, 2.24344807])
+    >>> sgd_dual_pi
+    array([[1.97377795e-02, 7.86706853e-02, 6.15682001e-03, 4.82586997e-09],
+           [4.19566963e-03, 4.42016865e-04, 2.02777272e-05, 3.68823708e-07],
+           [3.00379244e-03, 6.56562018e-02, 2.40462171e-02, 1.63579656e-10],
+           [4.28626062e-02, 5.60031599e-02, 6.13193826e-02, 9.67977735e-07],
+           [7.61972739e-02, 5.94609051e-03, 3.77886693e-02, 1.36046648e-03],
+           [4.44810042e-02, 1.89476742e-03, 1.73285847e-03, 3.51826036e-04],
+           [6.30118293e-02, 4.12398660e-02, 4.95148998e-02, 5.26247246e-06]])
 
     References
     ----------
diff --git a/ot/unbalanced.py b/ot/unbalanced.py
index bad12d6..50ec03c 100644
--- a/ot/unbalanced.py
+++ b/ot/unbalanced.py
@@ -6,6 +6,7 @@ Regularized Unbalanced OT
 # Author: Hicham Janati <hicham.janati@inria.fr>
 # License: MIT License
 
+from __future__ import division
 import warnings
 import numpy as np
 # from .utils import unif, dist
@@ -13,7 +14,7 @@ import numpy as np
 
 def sinkhorn_unbalanced(a, b, M, reg, alpha, method='sinkhorn', numItermax=1000,
                         stopThr=1e-9, verbose=False, log=False, **kwargs):
-    u"""
+    r"""
     Solve the unbalanced entropic regularization optimal transport problem and return the loss
 
     The function solves the following optimization problem:
@@ -75,7 +76,7 @@ def sinkhorn_unbalanced(a, b, M, reg, alpha, method='sinkhorn', numItermax=1000,
     >>> M=[[0., 1.], [1., 0.]]
     >>> ot.sinkhorn_unbalanced(a, b, M, 1, 1)
     array([[0.51122823, 0.18807035],
-          [0.18807035, 0.51122823]])
+           [0.18807035, 0.51122823]])
 
 
     References
@@ -122,7 +123,7 @@ def sinkhorn_unbalanced(a, b, M, reg, alpha, method='sinkhorn', numItermax=1000,
 def sinkhorn_unbalanced2(a, b, M, reg, alpha, method='sinkhorn',
                          numItermax=1000, stopThr=1e-9, verbose=False,
                          log=False, **kwargs):
-    u"""
+    r"""
     Solve the entropic regularization unbalanced optimal transport problem and return the loss
 
     The function solves the following optimization problem:
@@ -233,7 +234,7 @@ def sinkhorn_unbalanced2(a, b, M, reg, alpha, method='sinkhorn',
 
 def sinkhorn_knopp_unbalanced(a, b, M, reg, alpha, numItermax=1000,
                               stopThr=1e-9, verbose=False, log=False, **kwargs):
-    """
+    r"""
     Solve the entropic regularization unbalanced optimal transport problem and return the loss
 
     The function solves the following optimization problem:
@@ -287,12 +288,12 @@ def sinkhorn_knopp_unbalanced(a, b, M, reg, alpha, numItermax=1000,
     --------
 
     >>> import ot
-    >>> a=[.5, .15]
+    >>> a=[.5, .5]
     >>> b=[.5, .5]
     >>> M=[[0., 1.],[1., 0.]]
-    >>> ot.sinkhorn_knopp_unbalanced(a, b, M, 1., 1.)
-    array([[0.52761554, 0.22392482],
-           [0.10286295, 0.32257641]])
+    >>> ot.unbalanced.sinkhorn_knopp_unbalanced(a, b, M, 1., 1.)
+    array([[0.51122823, 0.18807035],
+           [0.18807035, 0.51122823]])
 
     References
     ----------
@@ -401,7 +402,7 @@ def sinkhorn_knopp_unbalanced(a, b, M, reg, alpha, numItermax=1000,
 
 def barycenter_unbalanced(A, M, reg, alpha, weights=None, numItermax=1000,
                           stopThr=1e-4, verbose=False, log=False):
-    """Compute the entropic regularized unbalanced wasserstein barycenter of distributions A
+    r"""Compute the entropic regularized unbalanced wasserstein barycenter of distributions A
 
      The function solves the following optimization problem:
 
diff --git a/ot/utils.py b/ot/utils.py
index efd1288..f21ceb9 100644
--- a/ot/utils.py
+++ b/ot/utils.py
@@ -285,9 +285,9 @@ class deprecated(object):
     The optional extra argument will be appended to the deprecation message
     and the docstring. Note: to use this with the default value for extra, put
     in an empty of parentheses:
-    >>> from ot.deprecation import deprecated
-    >>> @deprecated()
-    ... def some_function(): pass
+    >>> from ot.deprecation import deprecated  # doctest: +SKIP
+    >>> @deprecated()  # doctest: +SKIP
+    ... def some_function(): pass  # doctest: +SKIP
 
     Parameters
     ----------