Project CSP patterns to source

[SophieHerbst]

Hello! I am trying to implement Britta's and Jean-Remy's best of both world approach to project the CSP decoding patterns to source. I am able to do it on my data in my own code, but I am struggling to extract the patterns and coefficients in the pipeline code. On the one hand, I am a bit lost in the terminology between patterns, coefficients, and filters. On the other hand, when using Britta's code, I get an error message:

# Classifier
    csp = CSP(
        n_components=4,  # XXX revisit
        reg=0.1,  # XXX revisit
    )
    clf = make_pipeline(
        *preproc_steps,
        csp,
        LogReg(
            solver="liblinear",  # much faster than the default
            random_state=cfg.random_state,
            n_jobs=1,
        ),
    )

clf.fit(X, y)
        weights_csp = mne.decoding.get_coef(clf, 'patterns_', inverse_transform=True)

This estimator does not have a patterns_ attribute:
LogReg(n_jobs=1, random_state=42, solver='liblinear')

Can you help me understand how I need to change the decoding pipeline to make it work?

[britta-wstnr]

Hey @SophieHerbst, thanks for pinging me via email. I'd have to look into this a bit more - I'll put it on my agenda but please do ping me again if I don't move on it for a while, my to do list is a bit populated right now. One first thing to share: please note that the analyses for this paper were done quite a long time ago (around 2016, 2017 I would say). So this might also come down to changes in versions. But I will have a more detailed look!

[britta-wstnr]

PS: Cool that you want to add this! :)

[SophieHerbst]

PS: Cool that you want to add this! :)

if @larsoner and @hoechenberger agree of course. But to me this is the missing piece I have been looking for in a while to be able to interpret the CSP results.

[hoechenberger]

if @larsoner and @hoechenberger agree of course

👍👍👍👍👍

[larsoner]

Yes so far I have been projecting inverse-transformed CSP clf patterns_ to source space separately, the pipeline should do it for us. But given the pipeline runs all sensor analyses then all source analyses, we'll need a new source step whose job it is (just) to visualize the CSP patterns in source space

[SophieHerbst]

I think I found it. I needed to add a LinearModel wrapper in the clf

[SophieHerbst]

With the above code, I get an error in _05_decoding_csp.py, at csp.fit_transform() Do you have any idea on how to fix this?

A critical error occurred. The error message was: SVD did not converge

Aborting pipeline run. The traceback is:

  File "/home/sh254795/anaconda3/envs/mne/lib/python3.10/site-packages/mne_bids_pipeline/steps/sensor/_05_decoding_csp.py", line 354, in one_subject_decoding
    csp.fit_transform(X, y)
  File "/home/sh254795/anaconda3/envs/mne/lib/python3.10/site-packages/mne/decoding/csp.py", line 250, in fit_transform
    return super().fit_transform(X, y=y, **fit_params)
  File "/home/sh254795/anaconda3/envs/mne/lib/python3.10/site-packages/mne/decoding/mixin.py", line 33, in fit_transform
    return self.fit(X, y, **fit_params).transform(X)
  File "/home/sh254795/anaconda3/envs/mne/lib/python3.10/site-packages/mne/decoding/csp.py", line 197, in fit
    self.patterns_ = pinv(eigen_vectors)
  File "/home/sh254795/anaconda3/envs/mne/lib/python3.10/site-packages/mne/fixes.py", line 869, in pinv
    u, s, vh = np.linalg.svd(a, full_matrices=False)
  File "<__array_function__ internals>", line 180, in svd
  File "/home/sh254795/anaconda3/envs/mne/lib/python3.10/site-packages/numpy/linalg/linalg.py", line 1657, in svd
    u, s, vh = gufunc(a, signature=signature, extobj=extobj)
  File "/home/sh254795/anaconda3/envs/mne/lib/python3.10/site-packages/numpy/linalg/linalg.py", line 98, in _raise_linalgerror_svd_nonconvergence
    raise LinAlgError("SVD did not converge")

[larsoner]

I sometimes hit these and try to open issues on the relevant upstream trackers. In this case it's probably an OpenBLAS bug. You can try setting OMP_NUM_THREADS=1 in your env, sometimes it fixes it. In principle using mne.fixes._safe_svd could fix it but you'd have to be careful about broadcasting (it only works on arrays of shape (M, N) whereas np.linalg.svd works the last two dimensions of any array, i.e., those of shape (..., M, N).

Can you run the pipeline with --pdb, and do a np.save of a and print some info like:

(pdb) p np.savez("problem.npz", a)
(pdb) p gufunc
...
(pdb) p signature
...
(pdb) p extobj
...

and then also do (pdb) p np.show_config() and look to see where BLAS/LAPACK info exists, like:

>>> np.show_config()
Build Dependencies:
  blas:
    ...
    openblas configuration: OpenBLAS 0.3.27  USE64BITINT DYNAMIC_ARCH NO_AFFINITY

and upload the .npz and outputs here? Then I can try to replicate and escalate to the right place

[larsoner]

... there is also a chance that upgrading your OpenBLAS / NumPy could fix the issue. Or switching away from using MKL if you're using it

[larsoner]

... in other words this sounds a lot like https://github.com/OpenMathLib/OpenBLAS/issues/3044 . A couple of other useful outputs would be:

$ cat /sys/devices/cpu/caps/pmu_name
skylake
$ cat /proc/cpuinfo | grep 'model name' | uniq 
model name  : Intel(R) Core(TM) i7-7700K CPU @ 4.20GHz
$ mne sys_info | grep numpy
├☑ numpy             2.1.0.dev0+git20240514.35c4f37 (OpenBLAS 0.3.27 with 1 thread)

[SophieHerbst]

Thanks @larsoner The system specs are: cat /sys/devices/cpu/caps/pmu_name skylake cat /proc/cpuinfo | grep 'model name' | uniq model name : Intel(R) Xeon(R) Gold 5218 CPU @ 2.30GHz mne sys_info | grep numpy ├☑ numpy 1.23.5 (OpenBLAS 0.3.21 with 64 threads)

I updated numpy and am trying to run it again right now. ├☑ numpy 1.26.4 (OpenBLAS 0.3.23.dev with 64 threads)