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

    Computes :math:`s = \sum_t x_t` and :math:`C = X^\top X` while exploiting
    zero or constant columns in the data matrix.

    Parameters
    ----------
    X : ndarray (T, M)
        Data matrix
    remove_mean : bool
        True: remove column mean from the data, False: don't remove mean.
    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.
    weights: None or ndarray(T, )
        weights assigned to each trajectory point. If None, all data points have weight one.
        If ndarray, each data point is assigned a separate weight.
    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 : ndarray (M)
        sum
    C : ndarray (M, M)
        unnormalized covariance matrix

    """
    # Check consistency of inputs:
    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'
    # sparsify
    X0, mask_X, xconst = _sparsify(X, remove_mean=remove_mean, modify_data=modify_data,
                                   sparse_mode=sparse_mode, sparse_tol=sparse_tol)
    is_sparse = mask_X is not None
    # copy / convert
    # TODO: do we need to copy xconst?
    X0, xconst = _copy_convert(X0, const=xconst, remove_mean=remove_mean,
                               copy=is_sparse or (remove_mean and not modify_data))
    # sum / center
    w, sx, sx0_centered = _sum(X0, xmask=mask_X, xconst=xconst, symmetric=False, remove_mean=remove_mean,
                               weights=weights)
    if remove_mean:
        _center(X0, w, sx, mask=mask_X, const=xconst, inplace=True)  # fast in-place centering
    # TODO: we could make a second const check here. If after summation not enough zeros have appeared in the
    # TODO: consts, we switch back to dense treatment here.
    # compute covariance matrix
    if column_selection is not None:
        if is_sparse:
            Xk = X[:, column_selection]
            mask_Xk = mask_X[column_selection]
            X0k = Xk[:, mask_Xk]
            xksum = sx0_centered[column_selection]
            xkconst = Xk[0, ~mask_Xk]
            X0k, xkconst = _copy_convert(X0k, const=xkconst, remove_mean=remove_mean,
                                         copy=True)
            C = _M2(X0, X0k, mask_X=mask_X, mask_Y=mask_Xk, xsum=sx0_centered, xconst=xconst, ysum=xksum, yconst=xkconst,
                    weights=weights)
        else:
            X0k = X0[:, column_selection]
            C = _M2(X0, X0k, mask_X=mask_X, mask_Y=mask_X, xsum=sx0_centered, xconst=xconst,
                    ysum=sx0_centered[column_selection], yconst=xconst, weights=weights)
    else:
        C = _M2(X0, X0, mask_X=mask_X, mask_Y=mask_X, xsum=sx0_centered, xconst=xconst, ysum=sx0_centered, yconst=xconst,
                weights=weights, diag_only=diag_only)
    return w, sx, C
r""" Computes the first two unnormalized moments of X

    Computes :math:`s = \sum_t x_t` and :math:`C = X^\top X` while exploiting
    zero or constant columns in the data matrix.

    Parameters
    ----------
    X : ndarray (T, M)
        Data matrix
    remove_mean : bool
        True: remove column mean from the data, False: don't remove mean.
    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.
    weights: None or ndarray(T, )
        weights assigned to each trajectory point. If None, all data points have weight one.
        If ndarray, each data point is assigned a separate weight.
    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 : ndarray (M)
        sum
    C : ndarray (M, M)
        unnormalized covariance matrix