def OT_mapping_linear(xs, xt, reg=1e-6, ws=None,
                      wt=None, bias=True, log=False):
    """ return OT linear operator between samples

    The function estimates the optimal linear operator that aligns the two
    empirical distributions. This is equivalent to estimating the closed
    form mapping between two Gaussian distributions :math:`N(\mu_s,\Sigma_s)`
    and :math:`N(\mu_t,\Sigma_t)` as proposed in [14] and discussed in remark 2.29 in [15].

    The linear operator from source to target :math:`M`

    .. math::
        M(x)=Ax+b

    where :

    .. math::
        A=\Sigma_s^{-1/2}(\Sigma_s^{1/2}\Sigma_t\Sigma_s^{1/2})^{1/2}
        \Sigma_s^{-1/2}
    .. math::
        b=\mu_t-A\mu_s

    Parameters
    ----------
    xs : np.ndarray (ns,d)
        samples in the source domain
    xt : np.ndarray (nt,d)
        samples in the target domain
    reg : float,optional
        regularization added to the diagonals of convariances (>0)
    ws : np.ndarray (ns,1), optional
        weights for the source samples
    wt : np.ndarray (ns,1), optional
        weights for the target samples
    bias: boolean, optional
        estimate bias b else b=0 (default:True)
    log : bool, optional
        record log if True


    Returns
    -------
    A : (d x d) ndarray
        Linear operator
    b : (1 x d) ndarray
        bias
    log : dict
        log dictionary return only if log==True in parameters


    References
    ----------

    .. [14] Knott, M. and Smith, C. S. "On the optimal mapping of
        distributions", Journal of Optimization Theory and Applications
        Vol 43, 1984

    .. [15]  Peyré, G., & Cuturi, M. (2017). "Computational Optimal
        Transport", 2018.


    """

    d = xs.shape[1]

    if bias:
        mxs = xs.mean(0, keepdims=True)
        mxt = xt.mean(0, keepdims=True)

        xs = xs - mxs
        xt = xt - mxt
    else:
        mxs = np.zeros((1, d))
        mxt = np.zeros((1, d))

    if ws is None:
        ws = np.ones((xs.shape[0], 1)) / xs.shape[0]

    if wt is None:
        wt = np.ones((xt.shape[0], 1)) / xt.shape[0]

    Cs = (xs * ws).T.dot(xs) / ws.sum() + reg * np.eye(d)
    Ct = (xt * wt).T.dot(xt) / wt.sum() + reg * np.eye(d)

    Cs12 = linalg.sqrtm(Cs)
    Cs_12 = linalg.inv(Cs12)

    M0 = linalg.sqrtm(Cs12.dot(Ct.dot(Cs12)))

    A = Cs_12.dot(M0.dot(Cs_12))

    b = mxt - mxs.dot(A)

    if log:
        log = {}
        log['Cs'] = Cs
        log['Ct'] = Ct
        log['Cs12'] = Cs12
        log['Cs_12'] = Cs_12
        return A, b, log
    else:
        return A, b
return OT linear operator between samples

    The function estimates the optimal linear operator that aligns the two
    empirical distributions. This is equivalent to estimating the closed
    form mapping between two Gaussian distributions :math:`N(\mu_s,\Sigma_s)`
    and :math:`N(\mu_t,\Sigma_t)` as proposed in [14] and discussed in remark 2.29 in [15].

    The linear operator from source to target :math:`M`

    .. math::
        M(x)=Ax+b

    where :

    .. math::
        A=\Sigma_s^{-1/2}(\Sigma_s^{1/2}\Sigma_t\Sigma_s^{1/2})^{1/2}
        \Sigma_s^{-1/2}
    .. math::
        b=\mu_t-A\mu_s

    Parameters
    ----------
    xs : np.ndarray (ns,d)
        samples in the source domain
    xt : np.ndarray (nt,d)
        samples in the target domain
    reg : float,optional
        regularization added to the diagonals of convariances (>0)
    ws : np.ndarray (ns,1), optional
        weights for the source samples
    wt : np.ndarray (ns,1), optional
        weights for the target samples
    bias: boolean, optional
        estimate bias b else b=0 (default:True)
    log : bool, optional
        record log if True


    Returns
    -------
    A : (d x d) ndarray
        Linear operator
    b : (1 x d) ndarray
        bias
    log : dict
        log dictionary return only if log==True in parameters


    References
    ----------

    .. [14] Knott, M. and Smith, C. S. "On the optimal mapping of
        distributions", Journal of Optimization Theory and Applications
        Vol 43, 1984

    .. [15]  Peyré, G., & Cuturi, M. (2017). "Computational Optimal
        Transport", 2018.