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]

Typically this is all I use. Why also take the csp.patterns_? csp.fit(X, y) results anywhere else so it's a bit weird to me to add them here. It seems like we should just use these weights_csp. And maybe we should call them clf_patterns_ because they're really the patterns inverse transformed from the CSP all the way back through the other steps (e.g., PCA, sensor scaling)?

[SophieHerbst]

that is exactly where I am lost in the terminology..

my approach was to save the weights in this step to pick them up in a separate source projection step later

[SophieHerbst]

I wanted to save both, because in @britta-wstnr's code they are combined during beamforming:

stc_csp = beamform_components(weights_csp, sensor_pattern_csp, spat_filter,
                                fwd, multipliers=multiplier)

[britta-wstnr]

The idea here is to combine the coefficients and the CSP components. The (Haufe et al.) patterns are not easily beamform-able as they do not represent classical sink-source patterns anymore but refer to the decoding classes. Using the CSP components is a work-around for that.

[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 OpenMathLib/OpenBLAS#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)