<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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