How to use the pyriemann.utils.mean.mean_covariance function in pyriemann

# window length for calculating thresholds
    N = int(np.round(win_len * sfreq))

    if method == 'euclid':
        U = np.zeros((blocksize, nc, nc))
        for k in range(blocksize):
            rangevect = np.minimum(ns - 1, np.arange(k, ns + k, blocksize))
            x = X[:, rangevect]
            U[k, ...] = x @ x.T
        Uavg = block_geometric_median(U.reshape((-1, nc * nc)) / blocksize, 2)
        Uavg = Uavg.reshape((nc, nc))
    elif method == 'riemann':
        blocksize = int(ns // blocksize)
        U = block_covariance(X, window=blocksize, overlap=win_overlap)
        Uavg = pyriemann.utils.mean.mean_covariance(U, metric='riemann')

    # get the mixing matrix M
    M = linalg.sqrtm(np.real(Uavg))
    D, Vtmp = linalg.eig(M)
    # D, Vtmp = nonlinear_eigenspace(M, nc)
    V = Vtmp[:, np.argsort(D)]

    # get the threshold matrix T
    x = np.abs(np.dot(V, X))
    offsets = np.int_(np.arange(0, ns - N, np.round(N * (1 - win_overlap))))

    mu = np.zeros(nc)
    sig = np.zeros(nc)
    for ichan in range(nc):
        rms = x[ichan, :] ** 2
        Y = []

def test_mean_covariance_riemann():
    """Test mean_covariance for riemannian metric"""
    covmats, diags, A = generate_cov(100, 3)
    C = mean_covariance(covmats, metric='riemann')
    Ctrue = mean_riemann(covmats)
    assert_array_equal(C, Ctrue)

def __init_transform(self, X):
        """Init tr"""
        self.mdm = MDM(metric=self.metric, n_jobs=self.n_jobs)
        if self.mode == 'ftest':
            self.global_mean = mean_covariance(X, metric=self.mdm.metric_mean)
        elif self.mode == 'pairwise':
            X = pairwise_distance(X, metric=self.mdm.metric_dist)**2
        return X

sample_weight : None | ndarray shape (n_trials, 1)
            the weights of each sample. if None, each sample is treated with
            equal weights.

        Returns
        -------
        self : MDM instance
            The MDM instance.
        """
        self.classes_ = numpy.unique(y)

        if sample_weight is None:
            sample_weight = numpy.ones(X.shape[0])

        if self.n_jobs == 1:
            self.covmeans_ = [mean_covariance(X[y == l], metric=self.metric_mean,
                                    sample_weight=sample_weight[y == l])
                                        for l in self.classes_]
            """
            for l in self.classes_:
                self.covmeans_.append(
                    mean_covariance(X[y == l], metric=self.metric_mean,
                                    sample_weight=sample_weight[y == l]))
            """
        else:
            self.covmeans_ = Parallel(n_jobs=self.n_jobs)(
                delayed(mean_covariance)(X[y == l], metric=self.metric_mean,
                                         sample_weight=sample_weight[y == l])
                for l in self.classes_)

        return self

clean : array, shape=([n_trials, ]n_channels, n_samples)
        Clean data.
    state : 3-tuple
        Output ASR parameters.

    """
    M, T, R = state.values()
    [nc, ns] = X.shape

    if cov is None:
        if detrend:
            X_filt = signal.detrend(X_filt, axis=1, type='constant')
        cov = np.cov(X_filt, bias=True)
    else:
        if cov.ndim == 3:
            cov = pyriemann.utils.mean.mean_covariance(
                cov, metric='riemann', sample_weight=sample_weight)

    maxdims = int(np.fix(0.66 * nc))  # constant TODO make param

    # do a PCA to find potential artifacts
    if method == 'riemann':
        D, Vtmp = nonlinear_eigenspace(cov, nc)  # TODO
    else:
        D, Vtmp = np.linalg.eig(cov)

    V = np.real(Vtmp[:, np.argsort(D)])
    D = np.real(D[np.argsort(D)])

    # determine which components to keep (variance below directional threshold
    # or not admissible for rejection)
    keep = (D < np.sum(np.dot(T, V)**2, axis=0))

if sample_weight is None:
            sample_weight = numpy.ones(X.shape[0])

        if self.n_jobs == 1:
            self.covmeans_ = [mean_covariance(X[y == l], metric=self.metric_mean,
                                    sample_weight=sample_weight[y == l])
                                        for l in self.classes_]
            """
            for l in self.classes_:
                self.covmeans_.append(
                    mean_covariance(X[y == l], metric=self.metric_mean,
                                    sample_weight=sample_weight[y == l]))
            """
        else:
            self.covmeans_ = Parallel(n_jobs=self.n_jobs)(
                delayed(mean_covariance)(X[y == l], metric=self.metric_mean,
                                         sample_weight=sample_weight[y == l])
                for l in self.classes_)

        return self

def transform(self, X):
        """Tangent space projection.

        Parameters
        ----------
        X : ndarray, shape (n_trials, n_channels, n_channels)
            ndarray of SPD matrices.

        Returns
        -------
        ts : ndarray, shape (n_trials, n_ts)
            the tangent space projection of the matrices.
        """
        self._check_reference_points(X)
        if self.tsupdate:
            Cr = mean_covariance(X, metric=self.metric)
        else:
            Cr = self.reference_
        return tangent_space(X, Cr)

Parameters
        ----------
        X : ndarray, shape (n_trials, n_channels, n_channels)
            ndarray of SPD matrices.
        y : ndarray | None (default None)
            Not used, here for compatibility with sklearn API.
        sample_weight : ndarray | None (default None)
            weight of each sample.

        Returns
        -------
        self : TangentSpace instance
            The TangentSpace instance.
        """
        # compute mean covariance
        self.reference_ = mean_covariance(X, metric=self.metric,
                                          sample_weight=sample_weight)
        return self