This PR implements the discussion started in Issue #138

It consists of:

• A new sampling feature containing all the code relevant for sampling from Riemannian Gaussian distributions as defined in arXiv:1507.01760. I chose to give this rather general name because we could one day end up implementing other kinds of pdfs in the SPD manifold, like mixture of Gaussians or Riemannian Laplacian distributions. I am, however, open to suggestions!
• Two examples illustrating why this feature is useful. They reproduce the figures presented in Issue #138

It took me longer than initially expected because I ended up implementing a MCMC procedure for sampling the Riemannian Gaussian distribution. At least this ensures that we won't depend on other rather heavy packages like pyMC3 and pyro.

@sylvchev had mentioned the idea of also implementing the Bayesian classifiers from arXiv:1507.01760. This could indeed be useful for the package, however, I think it goes out of the scope of this PR. Furthermore, it does not apply directly the code that I'm including here... it depends mostly of an EM algorithm that I haven't touched yet. I suggest we keep this for (near) future work :)

Cheers, Pedro

Hi,

I'm trying to play with your code.

I was wondering whether it would be in principle possible to project a single time point (Nchan x 1) on the tangent space. I was thinking of fitting TangentSpace on the complete epochs, but then assess a series of svm on each time point separately in the tangent space. Although I'm not sure this makes sense at all..

Thanks!

Hi, everyone.

I would like to propose that we think of a logo for our dearest package ! :)

I have no training whatsoever in design and/or logo creation, so I must say that I don't have any strong opinion of what we should or should not include in a logo. With that said, I have a first suggestion and would like to know what you guys think.

We could also include colors, such as in the two proposals here below.

Well, that's it. I would be happy to know whether:

1. The logo itself seems OK or if we should go onto something else
2. What you think of the colors? We can of course play with other combinations (with the help of https://coolors.co/)

Cheers, Pedro

Hi, everyone

This is my first try to integrate some transfer learning methods to pyRiemann. @sylvchev and Emmanuel Kalunga had already started contributing on this topic with their PR #177 on MDWM, but I think it would be nice for us to discuss more broadly the general lines of development for transfer learning in our beloved package.

@agramfort @qbarthelemy and myself had a Zoom meeting recently in which we sketched some aspects of how an API for transfer learning in pyRiemann should look like. I tried implementing these ideas and had to make some design choices. The most relevant things are:

• Our data is always described by a triplet (X, y, meta) where X has the SPD matrices from all domains and y their corresponding labels. The meta structure is a pandas dataframe with two columns: domain for indicating to which domain each point belongs to and target indicating for each sample whether it is part of the Target domain or not.
• I had to create a TLSplitter object capable of splitting the data into training/validation partition in cross-validation. We can not use simply the ones from sklearn because we need the information from meta
• The way I see things, transfer learning procedures should be split into two parts: first we transform the data points and then we run classification. With this in mind, I created two examples of transformers DCT and RCT and a new object called TLPipeline that allows us to run the transfer learning procedure with a given classifier of choice.

I am, of course, open to suggestions and remarks and would like to know what you guys think of this first proposal.

Cheers, Pedro

I notice things are kinda stale here on the maintenance side. Notably, there are a few open PRs that look ready to merge (or at least there hasn't been any indication to the contrary).

I understand/assume that @alexandrebarachant is rather busy and only sporadically active in this repo, and I wouldn't ask him to be a more active maintainer (his time is probably better spent on other things than maintenance) so I'd like to suggest two possible solutions:

1. Move the repo to the NeuroTechX org, where there's a community that could actively and happily support it.
2. Add more maintainers to this repo.
   • From looking at the contributor stats, it looks like there aren't really any clear candidates. But @sylvchev has several open PRs and seems to have a good understanding of the subject matter.

The reason for asking this is that there are more issues here I'd like to open PRs to fix, but I won't put in the effort if I'm unlikely to get a response/merge anytime soon. An example would be to move from Travis (which has recently dropped their free-forever CI for open source) to GitHub Actions (see my changes in https://github.com/NeuroTechX/eeg-notebooks/pull/24 for an example).

There are probably more issues at hand, like how to deal with the PyPI upload rights, but I'd just like to throw these two options out and see what people think about these options.

Hi Alexandre,

I was trying the cross-validation classification example code from the homepage on some toy EEG data (64 channels), but I keep running into the following error:

I tracked the error down to the mean_riemann method . logm(tmp) (line 61) is filled with NaNs.

My input data is in the following format (Ntrials, Nchannels, Nsamples) :

Let me know if you need more information.

Thanks!

This PR adds an example to compare several Riemannian classifiers on low-dimensional synthetic datasets, adapted to SPD matrices from https://scikit-learn.org/stable/auto_examples/classification/plot_classifier_comparison.html

@gabelstein

Hello,

We implemented a SKLearn wrapper around the Qiskit library to work with EEG data using Riemannian Geometry.

Could this work be considered as inside the scope of pyRiemann?

This includes block diagonal covariance estimation. Useful for tasks where multiple band-pass filters are applied to the original data and covariance between frequencies do not hold information.

I'm actually using your implementation of XDawn in one of my studies. I would like to include a proper citation to the source (for now I mention the URL). Github supports generating proper DOI's though, which are easier to cite.

Making use of the new kernel module, the new support vector machine classification.

Should I extend some example code somewhere to compare e.g. MDM with the RSVC?

Hello again! I created a new function for frequency band selection on the manifold using the class distinctiveness measure.

This function finds the frequency band with the large class distinctiveness for training data from a broad frequency band that the user input. For instance, if the user inputs 5-35Hz, the function adjusts the band and returns 8-12Hz as the most class-distinctive frequency band. I added its example code using a motor imagery dataset, so please take a look.

Using an optimized frequency band is one of the important aspects for enhancing oscillatory activity classification, so hopefully, this new function is helpful for the user:-)

Description

A Hermitian positive definite (HPD) matrix $M$ is defined as: $M = A + iB$ where $A$ is a symmetric positive definite (SPD) matrix and $B$ is a skew-symmetric matrix.

Cross-spectral matrices are HPD matrices, where real parts are co-spectral matrices and imaginary parts are quadrature spectra capturing phase information.

Providing richer features, some works have shown the interest of classifying HPD matrices rather than SPD ones. 2017 - Dehgan - Classification in Riemannian space An application to sleep EEG.pdf 2019 - Xu - Feature extraction from the Hermitian manifold for Brain-Computer Interfaces.pdf

pyRiemann should be adapted to process such complex features. The first obvious change is to update pyriemann.utils.base._matrix_operator, adding .conj(): D = (eigvecs * eigvals) @ np.swapaxes(eigvecs.conj(), -2, -1)

@plcrodrigues @sylvchev @ygerf

Example

"""HPD matrices classification by MDM"""

import numpy as np
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.model_selection import cross_val_score, KFold

from pyriemann.classification import MDM
from pyriemann.utils.covariance import cross_spectrum


###############################################################################


class CrossSpectra(BaseEstimator, TransformerMixin):
    """Estimation of cross-spectral matrices.

    Complex cross-spectral matrices are HPD matrices estimated as the spectrum
    covariance in the frequency domain. It returns a 4-d array with a
    cross-spectral matrix for each input and in each frequency bin of the
    Fourier transform.

    Parameters
    ----------
    window : int, default=128
        The length of the FFT window used for spectral estimation.
    overlap : float, default=0.75
        The percentage of overlap between window.
    fmin : float | None, default=None
        The minimal frequency to be returned.
    fmax : float | None, default=None
        The maximal frequency to be returned.
    fs : float | None, default=None
        The sampling frequency of the signal.

    Attributes
    ----------
    freqs_ : ndarray, shape (n_freqs,)
        If transformed, the frequencies associated to cospectra.
        None if ``fs`` is None.

    See Also
    --------
    Covariances
    Coherences

    References
    ----------
    .. [1] https://en.wikipedia.org/wiki/Cross-spectrum
    """

    def __init__(self, window=128, overlap=0.75, fmin=None, fmax=None,
                 fs=None):
        """Init."""
        self.window = _nextpow2(window)
        self.overlap = overlap
        self.fmin = fmin
        self.fmax = fmax
        self.fs = fs

    def fit(self, X, y=None):
        """Fit.

        Do nothing. For compatibility purpose.

        Parameters
        ----------
        X : ndarray, shape (n_matrices, n_channels, n_times)
            Multi-channel time-series.
        y : None
            Not used, here for compatibility with sklearn API.

        Returns
        -------
        self : CrossSpectra instance
            The CrossSpectra instance.
        """
        return self

    def transform(self, X):
        """Estimate cross-spectral matrices.

        Parameters
        ----------
        X : ndarray, shape (n_matrices, n_channels, n_times)
            Multi-channel time-series.

        Returns
        -------
        X_new : ndarray, shape (n_matrices, n_channels, n_channels, n_freqs)
            Cross-spectral matrices for each input and for each frequency bin.
        """
        X_new = []

        for i in range(len(X)):
            S, freqs = cross_spectrum(
                X[i],
                window=self.window,
                overlap=self.overlap,
                fmin=self.fmin,
                fmax=self.fmax,
                fs=self.fs)
            X_new.append(S)
        self.freqs_ = freqs

        return np.array(X_new)


###############################################################################
# MDM on HPD matrices

n_matrices, n_channels, n_times = 50, 4, 5000
data_eeg = np.random.randn(2 * n_matrices, n_channels, n_times)

CrossSp = CrossSpectra(window=128, overlap=0.5, fmin=1, fmax=32, fs=128)
data_spec = CrossSp.transform(data_eeg)

ftarget = 10
X = np.squeeze(data_spec[..., CrossSp.freqs_ == ftarget])
y = [0] * n_matrices + [1] * n_matrices


mdm = MDM(metric='riemann')
cv = KFold(n_splits=10, shuffle=True, random_state=42)
scores = cross_val_score(mdm, X, y, cv=cv, n_jobs=1)

API change

Class CospCovariances could be renamed CoSpectra and simply coded as:

class CoSpectra(CrossSpectra):
    """Estimation of co-spectral matrices.

    Co-spectral matrices are SPD matrices estimated as the real part of the
    spectrum covariance in the frequency domain. It returns a 4-d array with a
    o-spectral matrix for each input and in each frequency bin of the
    Fourier transform.

    Parameters
    ----------
    window : int, default=128
        The length of the FFT window used for spectral estimation.
    overlap : float, default=0.75
        The percentage of overlap between window.
    fmin : float | None, default=None
        The minimal frequency to be returned.
    fmax : float | None, default=None
        The maximal frequency to be returned.
    fs : float | None, default=None
        The sampling frequency of the signal.

    Attributes
    ----------
    freqs_ : ndarray, shape (n_freqs,)
        If transformed, the frequencies associated to cospectra.
        None if ``fs`` is None.

    See Also
    --------
    Covariances
    CrossSpectra
    """

    def transform(self, X):
        """Estimate co-spectral matrices.

        Parameters
        ----------
        X : ndarray, shape (n_matrices, n_channels, n_times)
            Multi-channel time-series.

        Returns
        -------
        X_new : ndarray, shape (n_matrices, n_channels, n_channels, n_freqs)
            Co-spectral matrices for each input and for each frequency bin.
        """
        X_new = super().transform(X)
        return X_new.real

Formatting style across the project is inconsistent. I would suggest to fix this issue by using Black. It works pretty well, and it is very easy to use.

Thank you very much for such an excellent job ! You have pointed out in your article the advantages of robustness based on Riemannian geometry and gave a sample example about MI. But when I tried to apply this method to all the data of the datasets, I found that the total accuracy rate was only close to 63%, there are many subjects' acc are 50%. Can you give some suggestion what can I do now ?

I recently saw that there have been some improvements to the coherence calculation (https://github.com/alexandrebarachant/pyRiemann/pull/79/commits/9a6dbeef03f61c7658c8902923d0f22d84b57c1c) and I thought that I would like to propose some further improvements to the covariance and cospectra calculation.

In brief, I have observed that using einsum for these two operations give a speed-up of one order of magnitude. Here is my current code, but let us discuss if and how we could add this to pyriemann:

def covariances(x: np.ndarray) -> np.ndarray:
    """Calculate covariances on epoched data

    Input dimensions must be epoch, samples, channels
    """
    n = x.shape[1]
    # TODO: watch for einsum, it does not promote!
    c = np.einsum('aji,ajk->aik', x, x) / (n - 1)
    return c

def cross_spectrum(x: np.ndarray,
                   nperseg=None, noverlap=None, *,
                   fs: float = 1,
                   detrend: Optional[str] = 'constant',
                   window: str = 'boxcar',
                   #return_onesided: bool = True,
                   ) -> Tuple[np.ndarray, np.ndarray]:

    # x should be of shape (epoch, channel, sample)
    if x.ndim == 1:
        # when x is 1D, assume that is just samples; add a single channel
        x = x[np.newaxis, np.newaxis, :]
    elif x.ndim == 2:
        # when x is 2D, no manipulation is needed
        pass
    else:
        raise ValueError('Expected 1D or 2D array')

    n_channels, n_samples = x.shape
    nperseg = nperseg or n_samples
    noverlap = noverlap or 0

    # create sliding epochs to get epoch, channel, sample
    x_epoched = epoch(x, nperseg, nperseg - noverlap, axis=1)
    n_epochs = x_epoched.shape[0]

    # Handle detrending and window functions
    w = get_window(window, nperseg)
    x_epoched = x_epoched * w[np.newaxis, np.newaxis, :]
    if detrend is not None:
        x_epoched = signaltools.detrend(x_epoched, type=detrend, axis=2)

    # Apply FFT on x last dimension, X will be (epoch, channel, freq)
    freqs = np.fft.fftfreq(nperseg, 1 / fs)
    X = np.fft.fft(x_epoched)  # FFT over the last axis (samples)

    # Do a Einstein sum that will be equivalent to the following commented code:
    #
    # ## Verbose implementation ##
    # Apply x multiplied by its complex conjugate for each frequency
    # This gives dimensions epoch, channel, channel, frequency
    # cxx = np.apply_along_axis(_xxh, 1, X)
    #
    # Reorder the axis to epoch, frequency, channel, channel
    # cxx = np.rollaxis(cxx, 3, start=1)
    #
    # Average over epochs, eliminating the epoch dimension to get frequency, channel, channel
    # cxx = cxx.mean(axis=0)
    # ## end of verbose implementation ##
    #
    #
    # Using np.einsum, we get 1 order of magnitude faster (10x faster!),
    # but it is more difficult to understand. First, let us understand what
    # np.einsum('i,j->ij', a, b) does:
    # It multiplies the first axis of the first input over the each
    # element of the second input along its first axis.
    # In other words: it multiplies each element in a with each element in b
    # In other words: it does a vector outer product
    #
    # For example:
    # >>> x = np.arange(0, 3); y = np.arange(10, 13)
    # >>> x, y
    # (array([0, 1, 2]), array([10, 11, 12]))
    # >>> np.einsum('i,j->ij', x, y)
    # array([[ 0,  0,  0],
    #        [10, 11, 12],
    #        [20, 22, 24]])
    #
    # Another example:
    # >>> x = np.arange(3) + 1j
    # >>> x
    # array([0.+1.j, 1.+1.j, 2.+1.j])
    # >>> np.einsum('i,j->ij', x, x.conj())
    # array([[1.+0.j, 1.+1.j, 1.+2.j],
    #        [1.-1.j, 2.+0.j, 3.+1.j],
    #        [1.-2.j, 3.-1.j, 5.+0.j]])
    #
    # Back to our case: This is what we want to do on
    # an (I x J) array with I channels and J frequencies:
    # for each frequency, calculate x @ x.T (vector outer product):
    # np.einsum('ik,jk->kij', X, X.conj())
    #
    # for the 3D case (that is, with epochs)
    # np.einsum('ijl,ikl->iljk', X, X.conj())
    # or, for a more verbose approach, say the indices represent the following:
    # e: epoch
    # c: channel
    # f: frequency
    # h: channel on the conjugate tranpose (this should be the same size as c)
    # Then, the operation can be rewritten as:
    # np.einsum('ecf,ehf->efch', x, x.conj())
    #
    # Finally, since the final step is to do a mean over the epochs, we can
    # sum the "e" axis (by dropping the "e" axis on the output) and divide by
    # the number of epochs:
    cxx = np.einsum('ecf,ehf->fch', X, X.conj()) / n_epochs

    return freqs, cxx