nbara/python-meegkit

🔧🧠 MEEGkit: MEG & EEG processing toolkit in Python

python-meegkit's Issues

ASR transform fails (unexpected zi shape)

Hello! I am trying to apply ASR to some EEG data recorded at our lab with EEGLAB. We are using MNE for the analyses and we wanted to incorporate ASR as a preprocessing step. I was currently running the asr example code:

# Get data array for ASR
raw_array = raw.get_data()

# Train ASR
asr = ASR(method='euclid')
train_idx = np.arange(0 * sfreq, 30 * sfreq, dtype=int) # 30s period
_, sample_mask = asr.fit(raw_array[:, train_idx])

# Apply a sliding window to the data
X = sliding_window(raw_array, window=int(sfreq), step=int(sfreq))
Y = np.zeros_like(X)
for i in range(X.shape[1]):
    Y[:, i, :] = asr.transform(X[:, i, :])

And I got the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-63-b37fc0c9f43f> in <module>
     11 Y = np.zeros_like(X)
     12 for i in range(X.shape[1]):
---> 13     Y[:, i, :] = asr.transform(X[:, i, :])

~/gitrepos/python-meegkit/meegkit/asr.py in transform(self, X, y, **kwargs)
    221         # Yulewalk-filtered data (optional).
    222         X_filt, self.zi_ = yulewalk_filter(
--> 223             X, sfreq=self.sfreq, ab=self.ab_, zi=self.zi_)
    224 
    225         if not self._fitted:

~/gitrepos/python-meegkit/meegkit/utils/asr.py in yulewalk_filter(X, sfreq, zi, ab, axis)
    317     if zi is None:
    318         zi = signal.lfilter_zi(B, A)
--> 319         out, zf = signal.lfilter(B, A, X, zi=zi[:, None] * X[:, 0], axis=axis)
    320     else:
    321         out, zf = signal.lfilter(B, A, X, zi=zi, axis=axis)

~/anaconda3/envs/mne/lib/python3.7/site-packages/scipy/signal/signaltools.py in lfilter(b, a, x, axis, zi)
   1907             return sigtools._linear_filter(b, a, x, axis)
   1908         else:
-> 1909             return sigtools._linear_filter(b, a, x, axis, zi)
   1910 
   1911 

ValueError: Unexpected shape for zi:  expected (32,8), found (8,32).

My understanding of the analysis is not very deep, but after looking around in scipy documentation and the repo's code I found that this 'zi' parameter is calculated automatically before being passed to the lfilter function. Being so, I am a bit confused by the fact that it has unexpected shape. Any ideas? Thanks!

Applying ZapLine to continuous data

Hi, thanks for your package.

I was wondering if you suggest segmenting continuous data before applying the ZapLine algorithm. I tried to use it on a rather large chunk of data and it failed while trying to allocate an array of 32.2 TiB.

One option I can think of would be to cut my data into arbitrary segments of a few seconds each, apply ZapLine, and then piece the segments back together. What do you think?

ZapLine improvements

Todo:

• Write ZapLine example
  --> see #52
• Make ZapLine faster (data from #47 was really slow to run) and less memory-hungry (see #56)
  --> see #57
• Implement : https://github.com/maciekszul/zapline_iterator (cc @maciekszul)
  --> see #52
• Check out Zapline-plus : https://www.biorxiv.org/content/10.1101/2021.10.18.464805v1 (cc @MariusKlug )

complexwarning when using asr?

i have 2 different calibration data, one of them causes the complexwarning to appear but the other doesnt. no error occurs regardless, but the calibration data that causes the complexwarning results in a weird looking reconstructed signal.

ValueError: Covariance matrices must be positive definite. Add regularization to avoid this error.

I run into this error for some of my EEG inputs when I run the meegkit.asr.asr_process():

ValueError: Covariance matrices must be positive definite. Add regularization to avoid this error.

Complete error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-199-a17002d5a7b8> in <module>
     52             am = 0
     53             for i in range(len(X_t_w)):
---> 54                 x_t_w = asrC.transform(X_t_w[i].T).T
     55                 # x_t_w = octave.asr_process(X_t_w[i].T, sfreq, asr_state).T
     56                 X_t = avaDM.LagGenerator(x_t_w, lags)

c:\~\python\python38\lib\site-packages\meegkit\asr.py in transform(self, X, y, **kwargs)
    242 
    243         # Clean data, using covariances weighted by sample_weight
--> 244         out, self.state_ = asr_process(X, X_filt, self.state_,
    245                                        cov=np.stack(self.cov_),
    246                                        method=self.method,

c:\~\python\python38\lib\site-packages\meegkit\asr.py in asr_process(X, X_filt, state, cov, detrend, method, sample_weight)
    577     else:
    578         if cov.ndim == 3:
--> 579             cov = pyriemann.utils.mean.mean_covariance(
    580                 cov, metric='riemann', sample_weight=sample_weight)
    581 

c:\~\python\python38\lib\site-packages\pyriemann\utils\mean.py in mean_covariance(covmats, metric, sample_weight, *args)
    330         C = metric(covmats, sample_weight=sample_weight, *args)
    331     else:
--> 332         C = mean_methods[metric](covmats, sample_weight=sample_weight, *args)
    333     return C
    334 

c:\~\python\python38\lib\site-packages\pyriemann\utils\mean.py in mean_riemann(covmats, tol, maxiter, init, sample_weight)
     59         for index in range(Nt):
     60             tmp = numpy.dot(numpy.dot(Cm12, covmats[index, :, :]), Cm12)
---> 61             J += sample_weight[index] * logm(tmp)
     62 
     63         crit = numpy.linalg.norm(J, ord='fro')

c:\~\python\python38\lib\site-packages\pyriemann\utils\base.py in logm(Ci)
     44 
     45     """
---> 46     return _matrix_operator(Ci, numpy.log)
     47 
     48 

c:\~\python\python38\lib\site-packages\pyriemann\utils\base.py in _matrix_operator(Ci, operator)
      8     """matrix equivalent of an operator."""
      9     if Ci.dtype.char in typecodes['AllFloat'] and not numpy.isfinite(Ci).all():
---> 10         raise ValueError("Covariance matrices must be positive definite. Add regularization to avoid this error.")
     11     eigvals, eigvects = scipy.linalg.eigh(Ci, check_finite=False)
     12     eigvals = numpy.diag(operator(eigvals))

ValueError: Covariance matrices must be positive definite. Add regularization to avoid this error.

Running the same code on Matlab using the same EEG does not give me the error.

Questions about TSPCA

Hi @nbara ,

Thanks for this cool package.

We have some MEG data recorded from a KIT system (160 channels of data + 3 reference channels) and would like to try using your TSPCA implementation for reducing environmental noise in the data. We perform our analyses in MNE so it's nice to be able to do this in python. @larsoner

I just had a couple of questions: 1) In the README document, it says that "TSCPA may give inaccurate results due to insufficient testing (contributions welcome!)" I was wondering if there had been any updates about this? 2) Are there any recommended parameters when using the meegkit.tspca.tsr() function (e.g. for the shifts and weights)?

Thank you very much!

Error allocating huge Memory

Greetings,

i am exploring different packages and methods for preprocessing of meg/eeg data and i am really thankful to have found this translation of matlab, because i only use python so far. Although i following the documentation, i run into an error that doesn't make much sense to me. At this point i am not sure if its because of a wrong format of the data i put into the dss_line() and i can only think about MNE causing this problem.

Background of the data:

• raw data has about 3gb size
• contains 30eeg, 102 meg, 204 grad channels
• read in the .fif file with python MNE

this is the error i get after running the function:

_numpy.core._exceptions.ArrayMemoryError: Unable to allocate 32.4 TiB for an array with shape (2111000, 2111000) and data type float64

32,4TiB seems to be so huge, that i can't really understand how the data is transformed to such a big array. i allready tried to use only the 30 eeg channels and allso changing parameters nkeep, blocksize or using the dss_line_iter() does not change the amount of memory shown in the error.

i appreciate some advice what i am missing and some information if may MNE is causing the problem.

Thank you !

Implement trial-masked robust detrending

https://www.sciencedirect.com/science/article/pii/S0165027021000157

It's almost already doable with the current code.

What is missing is a function that takes event times and the max and min of your wanted epoch (so what you want to mask) and output appropriate weights as a (n_channels * n_times) matrix.

pypi installer

Hello, I really like your implementation of Artifact Subspace Rejection and I would like to import it in my own package.
I however cannot use your package as a requirement in PyPI because it doest not allow git install.

Would you be kind enought to upload to PyPI the dist from the current version that seems to nicely build?

Projecting RESS components on new observatios of data

Hello!

We have been using this package's RESS implementation in our recent work with SSVEPs and it has been a great help, thank you for that!

I wanted to ask a question regarding the 'maps' object being output from RESS (refered to as 'mixing matrix' in the docstring) and what they have diferent from the pure eigenvectors coming out from the linear eigendecomposition. We would be interested to know whether it is feasible to compute RESS on a set of observations and the use the maps (or eigenvectors) to transform a new observation into the same space of those outputed by the function. For example:

• An array of shape (4000, 32, 40) is inputed to the RESS function keeping just the first component, thus obtaining an output of shape (4000, 1, 40)

• I have a new example of shape (4000, 32) that I want to convert to the RESS component space (4000, 1) using the same components computed from the first transformation. How should I do this?

Thank you very much.

Issue with cca_crossvalidate

AA, RR, SD = cca.cca_crossvalidate(raw_fin, env_fin, shifts=None, sfreq=1, surrogate=False, plot=False) outputs an empty AA matrix.

Error when using a different sfreq

I try to run the following class instance:

from meegkit import asr

asr.ASR(sfreq=64)

but it gives me the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-264-fafe23240154> in <module>
      1 from meegkit import asr
      2 
----> 3 asr.ASR(sfreq=64)

c:\~\python\python38\lib\site-packages\meegkit\asr.py in __init__(self, sfreq, cutoff, blocksize, win_len, win_overlap, max_dropout_fraction, min_clean_fraction, name, method, **kwargs)
    124             80.0, (sfreq / 2.0) - 1.0), sfreq / 2.0]) * 2.0 / sfreq
    125         M = np.array([3, 0.75, 0.33, 0.33, 1, 1, 3, 3])
--> 126         B, A = yulewalk(8, F, M)
    127         self.ab_ = (A, B)
    128         self.cov_ = []

c:\~\python\python38\lib\site-packages\meegkit\utils\asr.py in yulewalk(order, F, M)
    241             inc = (j - nb) / (ne - nb)
    242 
--> 243         Ht[0][nb:ne + 1] = np.array(inc * M[i + 1] + (1 - inc) * M[i])
    244         nb = ne + 1
    245 

ValueError: could not broadcast input array from shape (385) into shape (257)

Information about filter bank analysis (FBA) and task-related component analysis (TRCA)

Hi,

I am currently working in a human factors research team to implement an SSVEP-based BCI of high-speed visual evoked potentials.

We are already using your library for part of our pipeline such as RESS, but we would like to test the approach explained in this document as well:

"An SSVEP-based high-speed BCI using dry EEG electrodes"
https://www.nature.com/articles/s41598-018-32283-8

In short, they show how it is possible to overcome the CCA algorithm by using a filter bank analysis (FBA) as part of preprocessing and the activity-related component analysis (TRCA) as a classification algorithm for the target classes.

Does Python-meegkit include or do you intend to include the FBA or TRCA algorithms?

Thank you.

dimension error when detrending with higher order polynomial

I am trying to do a detrending of my raw data in two steps with a linear and a higher order polynomial detrending, like this:

y, _, _ = meegkit.detrend.detrend(x, 1)
y, _, _ = meegkit.detrend.detrend(y, 5)

where x is raw data of shape (289120, 64)
The first detrending works just fine but during the second one I get a ValueError in line 209 of the regress function, saying
that the array c could not be broadcast into the matrix b due to a mismatch in dimensions.
It seems the matrix b has the shape (n_channels, n_channels) while c has the shape (order, 1).

Maybe something has changed in the underlying scipy functions? I am using:

python 3.9.2
numpy 1.20.2
scipy 1.6.2

Thanks for the nice package btw:)
Cheers,
Ole

module 'pymanopt' has no attribute 'function'

I cloned current master aed858d and run the following command pip install .

I tested examples/example_asr.ipynb
asr = ASR(method='euclid') works just well
but when I do:
asr = ASR(method='riemann')
the asr_fit() function works well, apparently since ISSUE#18 solved by enabling non scm estimators was addressed then merged into main

But then during asr_transform()
i get the following error
module 'pymanopt' has no attribute 'function'

I removed the @pymanopt.function.Callable and tested it again. Without success
The backtrack seems to point to the trust_regions function of pymanopt. it seems that the parameter
MAX_INNER_ITER: "maximum inner iterations" is 0 or NaN

Thanks a lot :)

Edit: i followed the example_asr.ipynb

support python 3.8

LinAlgError: SVD did not converge in Linear Least Squares

I get the following error the first time I run the ASR instance definition:

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-17-f449c895e6b2> in <module>
----> 1 asrC = asr.ASR(sfreq=125)

c:\~\python\python38\lib\site-packages\meegkit\asr.py in __init__(self, sfreq, cutoff, blocksize, win_len, win_overlap, max_dropout_fraction, min_clean_fraction, name, method, **kwargs)
    124             80.0, (sfreq / 2.0) - 1.0), sfreq / 2.0]) * 2.0 / sfreq
    125         M = np.array([3, 0.75, 0.33, 0.33, 1, 1, 3, 3])
--> 126         B, A = yulewalk(8, F, M)
    127         self.ab_ = (A, B)
    128         self.cov_ = []

c:\~\python\python38\lib\site-packages\meegkit\utils\asr.py in yulewalk(order, F, M)
    260         (1 / 2, np.ones((1, int(n2 - 1))), np.zeros((1, int(n - n2)))),
    261         axis=None)
--> 262     A = polystab(denf(R, order))  # compute denominator
    263 
    264     # compute additive decomposition

c:\~\python\python38\lib\site-packages\meegkit\utils\asr.py in denf(R, na)
    521     Rhs = - R[na + 1:nr]
    522     A = np.concatenate(
--> 523         (1, linalg.lstsq(Rm, Rhs.T, rcond=None)[0].T), axis=None)
    524     return A

<__array_function__ internals> in lstsq(*args, **kwargs)

c:\~\python\python38\lib\site-packages\numpy\linalg\linalg.py in lstsq(a, b, rcond)
   2257         # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis
   2258         b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
-> 2259     x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
   2260     if m == 0:
   2261         x[...] = 0

c:\~\python\python38\lib\site-packages\numpy\linalg\linalg.py in _raise_linalgerror_lstsq(err, flag)
    107 
    108 def _raise_linalgerror_lstsq(err, flag):
--> 109     raise LinAlgError("SVD did not converge in Linear Least Squares")
    110 
    111 def get_linalg_error_extobj(callback):

LinAlgError: SVD did not converge in Linear Least Squares

To reproduce the error, please run:

from meegkit import asr
asrC = asr.ASR(sfreq=125)

And my numpy version is 1.18.5, running python==3.8.6.

Q: mcca inverse transform

Hi @nbara ,

Thanks for this nice package.

IIUC, mcca generate a matrix that project sensors onto canonical components.

Is there a reverse transform easily available e.g. project many subjects on canonical components, average, and project back on one?

I think this should be doable given that canonical components are orthogonal and thus invertible?

Thanks!

Impact of component number on topographies and time courses

I have a conceptual question - why does the number of components (the keep1 parameter) influence the outcome of the dss?
If I understood the approach correctly it is essentially a PCA that finds the components with maximal power relative to the bias filter. If that is the case, shouldn't the first components always look the same, no matter how many I keep in total? From my data I get very different results when changing keep1.

for example, these are the first 4 components after the dss1 with keep1=5:

and these are the first 4 components with keep1=15:

Is this behavior expected?

'float' object has no attribute 'ppf' when executing `normfit()`

Hi @gferraro2019, I tried running example_trca.py, and encounter this error:

Block 0: accuracy = 97.5, 	ITR = 301.3
Block 1: accuracy = 100.0, 	ITR = 319.3
Block 2: accuracy = 95.0, 	ITR = 286.3
Block 3: accuracy = 95.0, 	ITR = 286.3
Block 4: accuracy = 95.0, 	ITR = 286.3
Block 5: accuracy = 100.0, 	ITR = 319.3
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-9-7ca177a962e5> in <module>
    124 
    125 # Mean accuracy and ITR computation
--> 126 mu, _, muci, _ = normfit(accs, alpha_ci)
    127 print()
    128 print(f"Mean accuracy = {mu:.1f}%\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f}%)")  # noqa

<ipython-input-4-769a1af8404d> in normfit(data, ci)
     43     num = len(arr)
     44     avg, std_err = np.mean(arr), scipy.stats.sem(arr)
---> 45     h_int = std_err * t.ppf((1 + ci) / 2., num - 1)
     46     var = np.var(data, ddof=1)
     47     var_ci_upper = var * (num - 1) / (chi2.ppf((1 - ci) / 2, num - 1))

AttributeError: 'float' object has no attribute 'ppf'

Issue with Removing Line Noise in MEG Data Using dss_line_iter

Hi everyone,

I used the dss_line_iter function to remove line noise from EEG and MEG data. The process works well on EEG but little strange on MEG.

Here is the function I use for removing line noise from MNE raw objects:

def linenoise_clean(raw, fline):
    # raw : mne.raw
    data = raw.get_data().T 
    data = np.expand_dims(data, axis=2) # make mne.raw dss_line_iter readable 
    sfreq = raw.info['sfreq']

    out, _ = dss.dss_line_iter(data, fline, sfreq,nfft=400,n_iter_max=100)
    cleaned_raw = mne.io.RawArray(out.T.squeeze(), raw.info)# array2mne.raw
    
    return cleaned_raw

For EEG, this method effectively eliminates the line noise, as shown in the PSD plot of raw EEG(above) and cleaned EEG(below).

However, when applying the same function to MEG data, the results are quite unusual. The line noise in the 0-300 Hz range remains consistent with the original data, and there’s an emergence of new line noise above 300 Hz, which was not present in the original data.

Has anyone else encountered this issue when applying dss_line_iter to MEG data? Any insights or suggestions on how to tweak the parameters or use an alternative approach to address this anomaly would be greatly appreciated.

Thank you in advance for your help and suggestions!

IndexError with asr.fit() for a [100s-150s] window and sfreq=500 Hz

Hi,

Thank you for this great library! I am facing a weird problem while trying to use asr.fit(). Following your ASR example, below is the code I want to run:

start_window = 100 
end_window = 150
sfreq = 500
asr = ASR(sfreq=sfreq,method='riemann')
train_idx = np.arange(start_window * sfreq, end_window * sfreq, dtype=int)
_, sample_mask = asr.fit(EEG_filt_data[:, train_idx])

But I get the following error:

  File ".../EEG_preprocessing.py", line 367, in <module>
    _, sample_mask = asr.fit(EEG_filt_data[:, train_idx])

  File ".../anaconda3/envs/mne/lib/python3.8/site-packages/meegkit/asr.py", line 173, in fit
    clean, sample_mask = clean_windows(

  File ".../anaconda3/envs/mne/lib/python3.8/site-packages/meegkit/asr.py", line 393, in clean_windows
    clean = np.delete(X, sample_mask2remove, 1)

  File "<__array_function__ internals>", line 5, in delete

  File ".../anaconda3/envs/mne/lib/python3.8/site-packages/numpy/lib/function_base.py", line 4480, in delete
    keep[obj,] = False

IndexError: arrays used as indices must be of integer (or boolean) type

It seems that the problem occurs only in the following situation (I don't know why):

• clean window of50s length;
• start_window fixed at 100;
• sfreq fixed at 500.

Indeed when I change one variable while keeping fixed the others (giving for instance 50, 100 and 499) it works fine. Maybe it is due to some kind of multiple relationship?

Versions used:
numpy = 1.20.1
mne = 0.22.1

Thank you!

dss_line_iter error

Hi, thanks for your handy package!
firstly i am new to eeg processing and programming, hence sorry if its a stupid question i just cant get solved.

i use the dss_line_iter function to get rid of line noise artefact but get thrown following error:

UnboundLocalError: local variable 'mean_psd' referenced before assignment
in the following: File "/home/analysis/anaconda3/lib/python3.9/site-packages/meegkit/dss.py", line 316, in dss_line_iter
mean_psd_wospot = mean_psd.copy()

the dss_line function without iter works fine.

Thanks for help!

RESS output shape appears to be (samples, trials) instead of (samples, channels, trials)

Hello!

I've been trying the recent RESS implementation of the library, and I am a bit puzzled at the output. The documentation on the ress function specifies that the input data should be in (samples, channels, trials) shape, and so it's done in the example script.

After trying with my own data of shape (3500, 32, 40), I get an output of shape (3500, 40). I also inspected a bit the code of the example, and the output for the simulated (1000, 10, 50) shaped data is also a (1000, 50) shaped array. From my understanding of the MATLAB implementation, I assumed the output would be the same shape as the input, but with the RESS components applied. Maybe I am missing something or the output of the RESS function is intended to be projected to the data, but I did not find any clues looking around the code.

Thanks in advance!

Memory usage of `dss_line()`

Hi,

I came across this repository when I was looking for an implementation of ZAPline. Really nice that you took the effort to translate the original matlab package. I have played around with it on my data, and for most of it, it works really nicely.

However, one thing that worries me a bit, is the memory usage of the dss_line(). So far, I tried to run it on MEG data (300 channels x 1.8x10^6 samples), but there was no chance that my laptop (24GB) could offer enough memory for dss_line to finish. Only if I cropped the data considerably it worked.

Of course, eventually the scripts will run on an HPC, so memory shouldn't be a big problem there, still it is quite annoying that I can't run and check those scripts locally before exporting them to the hpc. So I was wondering whether there are some magic tricks that I can do to lower the memory footprint? Settings that I have missed or something? I tried to trace the memory consumption of the function (where I was sufficiently confident to not affect functionality), and could indeed reduce the usage here and there, but the peak is still higher than what my pc can handle.

If helpful, I can provide more information of any kind.

Question on dss_line() usage

Hi, its awesome that there is an implementation of the ZAPline algorithm in Python! I'm trying to use it, but I somehow fail to remove the noise components with it. I'm convinced that I must be doing something obvious wrong - any idea what it may be?

I'm using Elekta MEG data (322 MEG channels, after Signal Space Separation) that show an artifact at 60Hz stemming from a presentation display.
Here is its PSD:

Here's the code I am running:

>>> import mne
>>> from meegkit.dss import dss_line
# read in the data with mne-python
>>> raw_fname = Path(datadir) / f'sub-{subject}/meg' / f'sub-{subject}_task-memento_proc-sss_meg.fif'
>>> raw = mne.io.read_raw_fif(raw_fname)
>>> raw.load_data()
>>> data, artifact = dss_line(raw._data.T, fline=60, sfreq=1000)
Power of components removed by DSS: 0.00

Putting the data back into the MNE Raw Object and visualizing the psd plot shows an almost identical profile, with the 60Hz component being pretty much unaffected.

Any idea what I might be doing wrong here?
Thanks much in advance!

Dataset format, can you clarify it?

Hello, can you clarify about dataset format? for example_asr.py

My dataset has the next structure in txt or excel formats and put this data to numpy
ch1 ch2 ch3 ch4 ch5 ch6 ch7 ch8
1 1 1 1 1 1 1 1

But when I use example for ASR and put this data in numpy to Raw (line - 17 in example_asr.py)- I receive error

TRCA has Low Accuracy for Identifying 8 Targets in 16-23 Frequency Range?

I attempted to use the TRCA algorithm with SSVEP data corresponding to 8 targets, flickering at frequencies between 16 to 23 Hz with 1 Hz difference between each stimulus and without any phase offset. Is it appropriate to use meegkit's TRCA implementation in this situation where there is less targets and no phase offset? I used the same filter bank presented in https://nbara.github.io/python-meegkit/auto_examples/example_trca.html, but the accuracy on my data (which I have ensured is in the correct format and filtered between 15 to 50 Hz before passing it through the filter banks) stays below 25% when using TRCA. I have in total 128 labels for a total of 16 experiment blocks, where each time the TRCA training gets 15 blocks (120 labels) and needs to predict 8 labels.
When I use CCA (from sci-kit learn implementation) on the same data, I get 87 % accuracy on average per experiment block. The fact that CCA works well with data that is formatted exactly as the data I am passing into TRCA tells me that the problem is not in the data formatting. I would appreciate any help figuring out how to get TRCA working, if that is possible, in the given scenario.

Bug in detrending of epoched data

There seems to be a bug in the detrending function. When i use it on two raw (two dimensional) data, the result is as expected but when I use it on epoched data, the result is very strange - probably there is an issue with reshaping the matrix, see my example below:

import numpy
from matplotlib import pyplot as plt
from meegkit.detrend import detrend

time = numpy.arange(0, 10, 0.01)

# on raw data the detrending works as expected:
X = numpy.zeros([len(time), 3])
X[:, 0], X[:, 1], X[:, 2] = numpy.sin(time)+5, numpy.sin(time), numpy.sin(time)-5

X += numpy.random.randn(X.shape[0], X.shape[1])

fix, ax = plt.subplots(2)
for x in X.T:
    ax[0].plot(time, x)
X, _, _ = detrend(X, 1)
for x in X.T:
    ax[1].plot(time, x)

# but on epoched data the result looks strange:
X = numpy.zeros([len(time), 3, 100])
X[:, 0] = numpy.repeat(numpy.sin(time)[:, numpy.newaxis], 100, axis=1)+5
X[:, 1] = numpy.repeat(numpy.sin(time)[:, numpy.newaxis], 100, axis=1)
X[:, 2] = numpy.repeat(numpy.sin(time)[:, numpy.newaxis], 100, axis=1)-5
X += numpy.random.randn(X.shape[0], X.shape[1], X.shape[2])

fix, ax = plt.subplots(2)
for x in X.mean(axis=2).T:
    ax[0].plot(time, x)
X, _, _ = detrend(X, 1)
for x in X.mean(axis=2).T:
    ax[1].plot(time, x)