MRG, BUG: Fix bugs with envcorr #8658

larsoner · 2020-12-14T16:20:09Z

Closes #8649.

@SherazKhan @dengemann can you look? This suggests the code I adapted originally has some bugs. It doesn't drastically change results or anything, in general it looks a bit cleaner now.

doc/changes/latest.inc

larsoner · 2020-12-14T16:22:44Z

mne/connectivity/envelope.py

            else:
                label_data_orth = (label_data * data_conj_scaled).imag
+                np.abs(label_data_orth, out=label_data_orth)


First change: when orthogonalizing, use abs of the orthogonalized data rather than just (data * data_conj_scaled).imag). Not 100% sure this is right but conceptually the thing we're correlating with is abs(data) so it seems to make sense to take the abs here, too...

I think the original operated on np.log10(label_data_orth ** 2) so an abs was implicitly taken by the squaring operation.

I am not sure about the log10, but they do square and take the log-transform:
Citation from their supplementary information:
"The signal Y(t,f) is orthogonalized in the complex plane with respect to X(t,f). This results in a positive number |Y_orth_X(t,f)| which is then squared and log-transformed, and correlated with [|X(t,f)|." Source

log10 is the same as log up to a scale factor, and that scale factor goes away in the correlation so it doesn't really matter. But np.log is faster than np.log10 so why not!

... but that text does make it clear that it should result in a positive number even before the square and log transform, so I think this is correct now (and will continue to be on the next push where I change np.log10 to np.log)

I think abs was missing in the original code, because of the squaring in the next step.

larsoner · 2020-12-14T16:23:38Z

mne/connectivity/envelope.py

                label_data_orth -= np.mean(label_data_orth, axis=-1,
                                           keepdims=True)
                label_data_orth_std = np.linalg.norm(label_data_orth, axis=-1)
                label_data_orth_std[label_data_orth_std == 0] = 1
            # correlation is dot product divided by variances
-            corr[li] = np.dot(label_data_orth, data_mag_nomean[li])
-            corr[li] /= data_mag_std[li]
+            corr[li] = np.sum(label_data_orth * data_mag_nomean, axis=1)


Second change is what @Avoide talked about in the other issue, we need to correlate with the correct pairs. This (plus the log=True, absolute=False) makes our code match both what @Avoide suggested and what FieldTrip produces.

Second change is what @Avoide talked about in the other issue, we need to correlate with the correct pairs. This (plus the log=True, absolute=False) makes our code match both what @Avoide suggested and what FieldTrip produces.

We were using np.corrcoef but this is much faster

mne/connectivity/tests/test_envelope.py

agramfort · 2020-12-14T17:36:47Z

@SherazKhan @dengemann it's your code originally. Can you take a minute to have a lookt? thanks

dengemann · 2020-12-14T18:00:29Z

Thank you everyone for looking into this. This work is currently not loaded into my working memory. What I can say though is that I always had a feeling since our working session at Biomag2018 that the MNE implementation does something different than the prototype.
I'll have a look at the code now.

dengemann · 2020-12-14T18:01:54Z

mne/connectivity/envelope.py

            corr = (corr.T + corr) / 2.
+            corr.flat[::corr.shape[0] + 1] = 0.


Perhaps consider adding a comment for Python newcomers that this updates the diagonal efficiently.

dengemann · 2020-12-14T18:13:00Z

@SherazKhan can you also take a look?

larsoner · 2020-12-14T18:28:48Z

mne/connectivity/tests/test_envelope.py

    # This is the version of the code by Sheraz and Denis.
    # For this version (epochs, labels, time) must be -> (labels, time, epochs)
-    data = np.transpose(data, (1, 2, 0))
-    corr_mats = np.empty((data.shape[0], data.shape[0], data.shape[2]))
-    for index, label_data in enumerate(data):
-        label_data_orth = np.imag(label_data * (data.conj() / np.abs(data)))
-        label_data_orig = np.abs(label_data)
-        label_data_cont = np.transpose(
-            np.dstack((label_data_orig, np.transpose(label_data_orth,
-                                                     (1, 2, 0)))), (1, 2, 0))
-        corr_mats[index] = np.array([np.corrcoef(dat)
-                                     for dat in label_data_cont])[:, 0, 1:].T
-    corr_mats = np.transpose(corr_mats, (2, 0, 1))
-    corr = np.mean(np.array([(np.abs(corr_mat) + np.abs(corr_mat).T) / 2.
-                             for corr_mat in corr_mats]), axis=0)


I always had a feeling since our working session at Biomag2018 that the MNE implementation does something different than the prototype.

@dengemann IIRC the code here was the prototype that I used to code a more efficient version in MNE, so it should have replicated it properly. Unfortunately I had to change this prototype code to do something different in order to get it to follow FieldTrip and the paper definitions...

No problem, If we can come up with a more flexible set of options this is good for research on this method.
The comparisons in the previous discussion suggest that difference may be rather subtle.

I second, differences will be subtle and additional options will be worth comparing.

larsoner · 2020-12-14T19:31:16Z

mne/connectivity/tests/test_envelope.py

+                y_orth_x = (y * x_conj_scaled).imag
+                y_orth_x_mag = np.abs(y_orth_x)
+                # Estimate correlation
+                corr[ii, jj] += np.abs(np.corrcoef(x_mag, y_orth_x_mag)[0, 1])


@SherazKhan the actual critical change in the code below is from data_mag_nomean[li] which only used the lith label's data to correlate with all N labels label_data_orth series, and this code now uses data_mag_nomean, i.e., all N labels data to correlate with all N labels label_data_orth series.

It's probably easier to see if you look at the "simplified" version here. The old/master version of the code does the correlation with the jj original (magnitude) data that you could call y_mag:

corr[ii, jj] += np.abs(np.corrcoef(np.abs(epoch_data[jj]), y_orth_x)[0, 1])

Instead of this correlation, with (magnitude) of the iith data that is called x_mag here, which is now implemented in this PR, equivalently written as:

corr[ii, jj] += np.abs(np.corrcoef(np.abs(epoch_data[ii]), y_orth_x)[0, 1])

@larsoner indeed, this is more closer to the original Hipp version.

larsoner · 2020-12-14T19:53:28Z

Updated example looks reasonable:

https://24108-1301584-gh.circle-artifacts.com/0/dev/auto_examples/connectivity/plot_mne_inverse_envelope_correlation.html#compute-the-degree-and-plot-it

dengemann · 2020-12-14T19:58:25Z

Updated example looks reasonable:

I agree the change looks reasonable.
But I am still baffled that the connectivity matrix looks so noisy, apparently lacking structure.
In the early examples we always saw some structure.

larsoner · 2020-12-14T20:05:10Z

The one on master isn't much better -- I think it's mostly because now we .crop(0, 60) for speed. If we remove that it (takes much longer and) looks like this on this PR:

dengemann · 2020-12-14T20:06:48Z

@larsoner yeah it's not a blocker, it was the same concern on master.

SherazKhan · 2020-12-14T20:13:39Z

The example should look similar to the original Hipp NN 2012

SherazKhan · 2020-12-14T20:16:27Z

@dengemann and I was able to reproduce something similar using cam-can data

dengemann · 2020-12-14T20:22:05Z

@SherazKhan I think we should repeat those benchmarks soon :)

SherazKhan · 2020-12-14T20:38:02Z

I did some testing on my data, please see the following for comparison between various options:

larsoner · 2020-12-15T12:21:32Z

@SherazKhan to my eye the clearest result is imag abs, which is now the default actually (and it's better than original imag, which is the default in master and is no longer available). Do you agree?

SherazKhan · 2020-12-15T12:31:21Z

@SherazKhan to my eye the clearest result is imag abs, which is now the default actually (and it's better than original imag, which is the default in master and is no longer available). Do you agree?

@larsoner I agree.

larsoner · 2020-12-15T17:54:37Z

Thanks for looking, discussing, and testing @dengemann @SherazKhan !

* upstream/master: (42 commits) MRG, ENH: Add DICS bias tests (mne-tools#8610) MRG, BUG, ENH: Add window option (mne-tools#8662) BUG: Fix alpha for volumes (mne-tools#8663) MRG, BUG: Fix bugs with envcorr (mne-tools#8658) MRG, ENH: Progressbar for csd_morlet (mne-tools#8608) Render is necessary now (mne-tools#8657) VIZ: Fix head size (mne-tools#8651) MRG, MAINT: bump sphinxcontrib-bitex version (mne-tools#8653) MRG, MAINT: Improve server env (mne-tools#8656) BUG: Mayavi center (mne-tools#8644) VIZ, ENH: allow show/hide annotations by label (mne-tools#8624) Add regression test for EEGLAB data with a chanlocs struct (mne-tools#8647) FIX: scalar_bar (mne-tools#8643) MRG: Small fix to tutorial; rename plot_events ordinate label to "Event id"; improve some SSP docstrings (mne-tools#8612) MRG, ENH: make plot alignment use defaults for colors (mne-tools#8553) BUG: Fix passing of channel type (mne-tools#8638) FIX: fixed loop over norm PSF/CTF options (mne-tools#8636) MRG, BUG: Pass kwargs (mne-tools#8630) DOC: Clearer error message (mne-tools#8631) BUG: Fix number of labels (mne-tools#8629) ...

* upstream/master: (38 commits) MRG, ENH: Add DICS bias tests (mne-tools#8610) MRG, BUG, ENH: Add window option (mne-tools#8662) BUG: Fix alpha for volumes (mne-tools#8663) MRG, BUG: Fix bugs with envcorr (mne-tools#8658) MRG, ENH: Progressbar for csd_morlet (mne-tools#8608) Render is necessary now (mne-tools#8657) VIZ: Fix head size (mne-tools#8651) MRG, MAINT: bump sphinxcontrib-bitex version (mne-tools#8653) MRG, MAINT: Improve server env (mne-tools#8656) BUG: Mayavi center (mne-tools#8644) VIZ, ENH: allow show/hide annotations by label (mne-tools#8624) Add regression test for EEGLAB data with a chanlocs struct (mne-tools#8647) FIX: scalar_bar (mne-tools#8643) MRG: Small fix to tutorial; rename plot_events ordinate label to "Event id"; improve some SSP docstrings (mne-tools#8612) MRG, ENH: make plot alignment use defaults for colors (mne-tools#8553) BUG: Fix passing of channel type (mne-tools#8638) FIX: fixed loop over norm PSF/CTF options (mne-tools#8636) MRG, BUG: Pass kwargs (mne-tools#8630) DOC: Clearer error message (mne-tools#8631) BUG: Fix number of labels (mne-tools#8629) ...

BUG: Fix bugs with envcorr

58d2eb1

larsoner commented Dec 14, 2020

View reviewed changes

larsoner added this to the 0.22 milestone Dec 14, 2020

larsoner changed the title ~~BUG: Fix bugs with envcorr~~ MRG, BUG: Fix bugs with envcorr Dec 14, 2020

DOC: Name, faster log

3176e3f

dengemann reviewed Dec 14, 2020

View reviewed changes

larsoner commented Dec 14, 2020

View reviewed changes

FIX: Invalid

b695a4d

larsoner commented Dec 14, 2020

View reviewed changes

FIX: Missed one

7a9839c

larsoner merged commit d4f99de into mne-tools:master Dec 15, 2020

larsoner deleted the envcorr branch December 15, 2020 17:54

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

MRG, BUG: Fix bugs with envcorr #8658

MRG, BUG: Fix bugs with envcorr #8658

larsoner commented Dec 14, 2020

larsoner Dec 14, 2020

Avoide Dec 14, 2020

larsoner Dec 14, 2020

larsoner Dec 14, 2020

SherazKhan Dec 14, 2020

larsoner Dec 14, 2020

SherazKhan Dec 14, 2020 •

edited

agramfort commented Dec 14, 2020

dengemann commented Dec 14, 2020

dengemann Dec 14, 2020

dengemann commented Dec 14, 2020

larsoner Dec 14, 2020

dengemann Dec 14, 2020

SherazKhan Dec 14, 2020 •

edited

larsoner Dec 14, 2020

SherazKhan Dec 14, 2020

larsoner commented Dec 14, 2020

dengemann commented Dec 14, 2020

larsoner commented Dec 14, 2020

dengemann commented Dec 14, 2020

SherazKhan commented Dec 14, 2020

SherazKhan commented Dec 14, 2020 •

edited

dengemann commented Dec 14, 2020

SherazKhan commented Dec 14, 2020 •

edited

larsoner commented Dec 15, 2020

SherazKhan commented Dec 15, 2020

larsoner commented Dec 15, 2020

		corr = (corr.T + corr) / 2.
		corr.flat[::corr.shape[0] + 1] = 0.

MRG, BUG: Fix bugs with envcorr #8658

MRG, BUG: Fix bugs with envcorr #8658

Conversation

larsoner commented Dec 14, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

SherazKhan Dec 14, 2020 • edited

Choose a reason for hiding this comment

agramfort commented Dec 14, 2020

dengemann commented Dec 14, 2020

Choose a reason for hiding this comment

dengemann commented Dec 14, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

SherazKhan Dec 14, 2020 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

larsoner commented Dec 14, 2020

dengemann commented Dec 14, 2020

larsoner commented Dec 14, 2020

dengemann commented Dec 14, 2020

SherazKhan commented Dec 14, 2020

SherazKhan commented Dec 14, 2020 • edited

dengemann commented Dec 14, 2020

SherazKhan commented Dec 14, 2020 • edited

larsoner commented Dec 15, 2020

SherazKhan commented Dec 15, 2020

larsoner commented Dec 15, 2020

SherazKhan Dec 14, 2020 •

edited

SherazKhan Dec 14, 2020 •

edited

SherazKhan commented Dec 14, 2020 •

edited

SherazKhan commented Dec 14, 2020 •

edited