def _M2_symmetric(Xvar, Yvar, mask_X=None, mask_Y=None, xsum=0, xconst=0, ysum=0, yconst=0, weights=None,
                  column_selection=None, diag_only=False):
    """ symmetric second moment matrices. Decide if we need dense, sparse, const"""
    if mask_X is None and mask_Y is None:
        if column_selection is None:
            Xvark = Xvar
            Yvark = Yvar
        else:
            Xvark = Xvar[:, column_selection]
            Yvark = Yvar[:, column_selection]
        Cxxyy = _M2_dense(Xvar, Xvark, weights=weights, diag_only=diag_only) \
                + _M2_dense(Yvar, Yvark, weights=weights, diag_only=diag_only)
        Cxy = _M2_dense(Xvar, Yvark, weights=weights, diag_only=diag_only)
        Cyx = _M2_dense(Yvar, Xvark, weights=weights, diag_only=diag_only)
        Cxyyx = Cxy + Cyx
    else:
        # Check if one of the masks is not None, modify it and also adjust the constant columns:
        if mask_X is None:
            mask_X = np.ones(Xvar.shape[1], dtype=np.bool)
            xconst = np.ones(0, dtype=float)
        if mask_Y is None:
            mask_Y = np.ones(Yvar.shape[1], dtype=np.bool)
            yconst = np.ones(0, dtype=float)
        if _is_zero(xsum) and _is_zero(ysum) or _is_zero(xconst) and _is_zero(yconst):
            Cxxyy, Cxyyx = _M2_sparse_sym(Xvar, mask_X, Yvar, mask_Y, weights=weights, column_selection=column_selection)
        else:
            xvarsum = xsum[mask_X]  # to variable part
            yvarsum = ysum[mask_Y]  # to variable part
            if column_selection is None:
                Xvark = Xvar
                mask_Xk = mask_X
                xkvarsum = xvarsum
                xkconst = xconst
                Yvark = Yvar
                mask_Yk = mask_Y
                ykvarsum = yvarsum
                ykconst = yconst
            else:
                Xvark = Xvar[:, _filter_variable_indices(mask_X, column_selection)]
                mask_Xk = mask_X[column_selection]
                xksum = xsum[column_selection]
                xkvarsum = xksum[mask_Xk]
                xkconst = xconst[_filter_variable_indices(~mask_X, column_selection)]
                Yvark = Yvar[:, _filter_variable_indices(mask_Y, column_selection)]
                mask_Yk = mask_Y[column_selection]
                yksum = ysum[column_selection]
                ykvarsum = yksum[mask_Yk]
                ykconst = yconst[_filter_variable_indices(~mask_Y, column_selection)]
            Cxxyy = _M2_const(Xvar, mask_X, xvarsum, xconst, Xvark, mask_Xk, xkvarsum, xkconst, weights=weights) \
                    + _M2_const(Yvar, mask_Y, yvarsum, yconst, Yvark, mask_Yk, ykvarsum, ykconst, weights=weights)
            Cxy = _M2_const(Xvar, mask_X, xvarsum, xconst, Yvark, mask_Yk, ykvarsum, ykconst, weights=weights)
            Cyx = _M2_const(Yvar, mask_Y, yvarsum, yconst, Xvark, mask_Xk, xkvarsum, xkconst, weights=weights)
            Cxyyx = Cxy + Cyx
    return Cxxyy, Cxyyx
symmetric second moment matrices. Decide if we need dense, sparse, const