def moments_XXXY(X, Y, remove_mean=False, symmetrize=False, weights=None,
                 modify_data=False, sparse_mode='auto', sparse_tol=0.0,
                 column_selection=None, diag_only=False):
    """ Computes the first two unnormalized moments of X and Y

    If symmetrize is False, computes

    .. math:
        s_x  &=& \sum_t x_t
        s_y  &=& \sum_t y_t
        C_XX &=& X^\top X
        C_XY &=& X^\top Y

    If symmetrize is True, computes

    .. math:
        s_x = s_y &=& \frac{1}{2} \sum_t(x_t + y_t)
        C_XX      &=& \frac{1}{2} (X^\top X + Y^\top Y)
        C_XY      &=& \frac{1}{2} (X^\top Y + Y^\top X)

    while exploiting zero or constant columns in the data matrix.

    Parameters
    ----------
    X : ndarray (T, M)
        Data matrix
    Y : ndarray (T, N)
        Second data matrix
    remove_mean : bool
        True: remove column mean from the data, False: don't remove mean.
    symmetrize : bool
        Computes symmetrized means and moments (see above)
    weights : None or ndarray(T, )
        weights assigned to each trajectory point of X. If None, all data points have weight one.
        If ndarray, each data point is assigned a separate weight.
    time_lagged : bool,
        indicates that Y is a time-lagged version of X.
    modify_data : bool
        If remove_mean=True, the mean will be removed in the data matrix X,
        without creating an independent copy. This option is faster but might
        lead to surprises because your input array is changed.
    sparse_mode : str
        one of:
            * 'dense' : always use dense mode
            * 'sparse' : always use sparse mode if possible
            * 'auto' : automatic
    sparse_tol: float
        Threshold for considering column to be zero in order to save computing
        effort when the data is sparse or almost sparse.
        If max(abs(X[:, i])) < sparse_tol, then row i (and also column i if Y
        is not given) of the covariance matrix will be set to zero. If Y is
        given and max(abs(Y[:, i])) < sparse_tol, then column i of the
        covariance matrix will be set to zero.
    column_selection: ndarray(k, dtype=int) or None
        Indices of those columns that are to be computed. If None, all columns are computed.
    diag_only: bool
        If True, the computation is restricted to the diagonal entries (autocorrelations) only.

    Returns
    -------
    w : float
        statistical weight
    s_x : ndarray (M)
        x-sum
    s_y : ndarray (N)
        y-sum
    C_XX : ndarray (M, M)
        unnormalized covariance matrix of X
    C_XY : ndarray (M, N)
        unnormalized covariance matrix of XY

    """
    # Check consistency of inputs:
    if Y is not None:
        assert Y.shape[0] == X.shape[0], 'X and Y must have equal length.'
    if weights is not None:
        assert X.shape[0] == weights.shape[0], 'X and weights_x must have equal length'
    # diag_only is only implemented for dense mode
    if diag_only and sparse_mode is not 'dense':
        if sparse_mode is 'sparse':
            import warnings
            warnings.warn('Computing diagonal entries only is not implemented for sparse mode. Switching to dense mode.')
        sparse_mode = 'dense'
    if diag_only and X.shape[1] != Y.shape[1]:
        raise ValueError('Computing diagonal entries only does not make sense for rectangular covariance matrix.')
    # sparsify
    X0, mask_X, xconst, Y0, mask_Y, yconst = _sparsify_pair(X, Y, remove_mean=remove_mean, modify_data=modify_data,
                                                            symmetrize=symmetrize, sparse_mode=sparse_mode, sparse_tol=sparse_tol)
    is_sparse = mask_X is not None and mask_Y is not None
    # copy / convert
    copy = is_sparse or (remove_mean and not modify_data)
    X0, xconst = _copy_convert(X0, const=xconst, remove_mean=remove_mean, copy=copy)
    Y0, yconst = _copy_convert(Y0, const=yconst, remove_mean=remove_mean, copy=copy)
    # sum / center
    w, sx, sx_centered, sy, sy_centered = _sum(X0, xmask=mask_X, xconst=xconst, Y=Y0, ymask=mask_Y, yconst=yconst,
                                               symmetric=symmetrize, remove_mean=remove_mean, weights=weights)
    if remove_mean:
        _center(X0, w, sx, mask=mask_X, const=xconst, inplace=True)  # fast in-place centering
        _center(Y0, w, sy, mask=mask_Y, const=yconst, inplace=True)  # fast in-place centering

    if symmetrize:
        Cxx, Cxy = _M2_symmetric(X0, Y0, mask_X=mask_X, mask_Y=mask_Y,
                                 xsum=sx_centered, xconst=xconst, ysum=sy_centered, yconst=yconst, weights=weights,
                                 column_selection=column_selection, diag_only=diag_only)
    else:
        if column_selection is not None:
            if is_sparse:
                Xk = X[:, column_selection]
                mask_Xk = mask_X[column_selection]
                X0k = Xk[:, mask_Xk]
                xksum = sx_centered[column_selection]
                xkconst = Xk[0, ~mask_Xk]
                X0k, xkconst = _copy_convert(X0k, const=xkconst, remove_mean=remove_mean,
                                             copy=True)

                Yk = Y[:, column_selection]
                mask_Yk = mask_Y[column_selection]
                Y0k = Yk[:, mask_Yk]
                yksum = sy_centered[column_selection]
                ykconst = Yk[0, ~mask_Yk]
                Y0k, ykconst = _copy_convert(Y0k, const=ykconst, remove_mean=remove_mean,
                                             copy=True)

                Cxx = _M2(X0, X0k, mask_X=mask_X, mask_Y=mask_Xk, xsum=sx_centered, xconst=xconst, ysum=xksum, yconst=xkconst,
                        weights=weights)
                Cxy = _M2(X0, Y0k, mask_X=mask_X, mask_Y=mask_Yk, xsum=sx_centered, xconst=xconst, ysum=yksum, yconst=ykconst,
                        weights=weights)
            else:
                X0k = X0[:, column_selection]
                Y0k = Y0[:, column_selection]
                Cxx = _M2(X0, X0k, mask_X=mask_X, mask_Y=mask_X, xsum=sx_centered, xconst=xconst,
                        ysum=sx_centered[column_selection], yconst=xconst, weights=weights)
                Cxy = _M2(X0, Y0k, mask_X=mask_X, mask_Y=mask_Y, xsum=sx_centered, xconst=xconst,
                        ysum=sy_centered[column_selection], yconst=yconst, weights=weights)
        else:
            Cxx = _M2(X0, X0, mask_X=mask_X, mask_Y=mask_X, xsum=sx_centered, xconst=xconst, ysum=sx_centered, yconst=xconst,
                      weights=weights, diag_only=diag_only)
            Cxy = _M2(X0, Y0, mask_X=mask_X, mask_Y=mask_Y, xsum=sx_centered, xconst=xconst, ysum=sy_centered, yconst=yconst,
                      weights=weights, diag_only=diag_only)

    return w, sx, sy, Cxx, Cxy
Computes the first two unnormalized moments of X and Y

    If symmetrize is False, computes

    .. math:
        s_x  &=& \sum_t x_t
        s_y  &=& \sum_t y_t
        C_XX &=& X^\top X
        C_XY &=& X^\top Y

    If symmetrize is True, computes

    .. math:
        s_x = s_y &=& \frac{1}{2} \sum_t(x_t + y_t)
        C_XX      &=& \frac{1}{2} (X^\top X + Y^\top Y)
        C_XY      &=& \frac{1}{2} (X^\top Y + Y^\top X)

    while exploiting zero or constant columns in the data matrix.

    Parameters
    ----------
    X : ndarray (T, M)
        Data matrix
    Y : ndarray (T, N)
        Second data matrix
    remove_mean : bool
        True: remove column mean from the data, False: don't remove mean.
    symmetrize : bool
        Computes symmetrized means and moments (see above)
    weights : None or ndarray(T, )
        weights assigned to each trajectory point of X. If None, all data points have weight one.
        If ndarray, each data point is assigned a separate weight.
    time_lagged : bool,
        indicates that Y is a time-lagged version of X.
    modify_data : bool
        If remove_mean=True, the mean will be removed in the data matrix X,
        without creating an independent copy. This option is faster but might
        lead to surprises because your input array is changed.
    sparse_mode : str
        one of:
            * 'dense' : always use dense mode
            * 'sparse' : always use sparse mode if possible
            * 'auto' : automatic
    sparse_tol: float
        Threshold for considering column to be zero in order to save computing
        effort when the data is sparse or almost sparse.
        If max(abs(X[:, i])) < sparse_tol, then row i (and also column i if Y
        is not given) of the covariance matrix will be set to zero. If Y is
        given and max(abs(Y[:, i])) < sparse_tol, then column i of the
        covariance matrix will be set to zero.
    column_selection: ndarray(k, dtype=int) or None
        Indices of those columns that are to be computed. If None, all columns are computed.
    diag_only: bool
        If True, the computation is restricted to the diagonal entries (autocorrelations) only.

    Returns
    -------
    w : float
        statistical weight
    s_x : ndarray (M)
        x-sum
    s_y : ndarray (N)
        y-sum
    C_XX : ndarray (M, M)
        unnormalized covariance matrix of X
    C_XY : ndarray (M, N)
        unnormalized covariance matrix of XY