Add TFR vector visualization

[alexrockhill]

Porting over mne-tools/mne-python#11524

MWE:

import numpy as np
import mne
import mne_gui_addons
from mne.datasets import somato
from mne.time_frequency import csd_tfr

data_path = somato.data_path()
subject = '01'
task = 'somato'
raw_fname = (data_path / 'sub-{}'.format(subject) / 'meg' /
             'sub-{}_task-{}_meg.fif'.format(subject, task))

# crop to 5 minutes to save memory
raw = mne.io.read_raw_fif(raw_fname).crop(0, 300)
raw.load_data()

picks = mne.pick_types(raw.info, meg='mag', eog=True, exclude='bads')

# read epochs
events = mne.find_events(raw)
epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, picks=picks,
                    preload=True, decim=3)

# read forward operator and point to freesurfer subject directory
fname_fwd = (data_path / 'derivatives' / 'sub-{}'.format(subject) /
             'sub-{}_task-{}-fwd.fif'.format(subject, task))
subjects_dir = data_path / 'derivatives' / 'freesurfer' / 'subjects'

fwd = mne.read_forward_solution(fname_fwd)

rank = mne.compute_rank(epochs, tol=1e-6, tol_kind='relative')
active_win = (0.5, 1.5)
baseline_win = (-1, 0)

# compute cross-spectral density matrices
freqs = np.logspace(np.log10(12), np.log10(30), 9)

# time-frequency decomposition
epochs_tfr = mne.time_frequency.tfr_morlet(
    epochs, freqs=freqs, n_cycles=freqs / 2, return_itc=False,
    average=False, output='complex')
epochs_tfr.decimate(20)  # decimate for speed

csd = csd_tfr(epochs_tfr, tmin=-1, tmax=1.5)
baseline_csd = csd_tfr(epochs_tfr, tmin=baseline_win[0], tmax=baseline_win[1])
ers_csd = csd_tfr(epochs_tfr, tmin=active_win[0], tmax=active_win[1])

surface = subjects_dir / subject / 'bem' / 'inner_skull.surf'
vol_src = mne.setup_volume_source_space(
    subject=subject, subjects_dir=subjects_dir, surface=surface,
    pos=10, add_interpolator=False)  # just for speed!

conductivity = (0.3,)  # one layer for MEG
model = mne.make_bem_model(subject=subject, ico=3,  # just for speed
                           conductivity=conductivity,
                           subjects_dir=subjects_dir)
bem = mne.make_bem_solution(model)

trans = fwd['info']['mri_head_t']
vol_fwd = mne.make_forward_solution(
    raw.info, trans=trans, src=vol_src, bem=bem, meg=True, eeg=True,
    mindist=5.0, n_jobs=1, verbose=True)

# Compute source estimate using MNE solver
snr = 3.0
lambda2 = 1.0 / snr ** 2
method = 'MNE'  # use MNE method (could also be dSPM or sLORETA)

# make a different inverse operator for each frequency so as to properly
# whiten the sensor data
inverse_operator = list()
for freq_idx in range(epochs_tfr.freqs.size):
    # for each frequency, compute a separate covariance matrix
    baseline_cov = baseline_csd.get_data(index=freq_idx, as_cov=True)
    baseline_cov['data'] = baseline_cov['data'].real  # only normalize by real
    # then use that covariance matrix as normalization for the inverse
    # operator
    inverse_operator.append(mne.minimum_norm.make_inverse_operator(
        epochs.info, vol_fwd, baseline_cov))

# finally, compute the stcs for each epoch and frequency
stcs = mne.minimum_norm.apply_inverse_tfr_epochs(
    epochs_tfr, inverse_operator, lambda2, method=method,
    pick_ori='vector')

viewer = mne_gui_addons.view_vol_stc(
    stcs, subject=subject, subjects_dir=subjects_dir,
    src=vol_src, inst=epochs_tfr)
viewer.go_to_extreme()  # show the maximum intensity source vertex
viewer.set_cmap(vmin=0.25, vmid=0.8)
viewer.set_3d_view(azimuth=40, elevation=35, distance=300)

[alexrockhill]

@larsoner looks like the examples don't quite build properly for 3D rendering. We could just remove them since they are redundant with mne-python or maybe not too much effort to fix: https://output.circle-artifacts.com/output/job/f775c93f-fe1b-47c9-a874-d1278d6f6097/artifacts/0/dev/auto_examples/evoked_ers_source_power.html#sphx-glr-auto-examples-evoked-ers-source-power-py

Also, I thought Black was supposed to fix your style for you, am I supposed to do something about the failed style test, sorry I'm not as familiar with that.

[larsoner]

For now you need to black . or so. We could add this as a pre-commit hook, or add a CI to do it

For the black rendering can you try adding this line

https://github.com/mne-tools/mne-python/blob/686857c55288ed91cc13cdc43bf100ab1f6c475d/.circleci/config.yml#L103

[alexrockhill]

Hmmm black . says all local files left unchanged but style still fails

[larsoner]

I'm guessing you're using a version older than 22.3.0, locally I have that and it changed some stuff so I pushed (along with a tiny commit to improve pyproject.toml completeness).

[alexrockhill]

Ah I think it was black installed as the debian and not the pipy distribution, just to note

[alexrockhill]

Also I don't have write access FYI

[alexrockhill]

@larsoner , have you solved this CI issue before?

[larsoner]

Yes PySide6 v6.5 broke some stuff, pushed a commit to hopefully fix it

[alexrockhill]

Hey all green, nice! I'm not quite sure what {0} argument is being passed to the CI so glad you knew!

[larsoner]

I think this is equivalent to data.real

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> data = rng.normal(size=(1000,)) + 1j * rng.normal(size=(1000,))
>>> np.allclose(data.real, np.abs(data) * np.cos(np.angle(data)))
True

[alexrockhill]

Interesting, so that was just plotting the real component, that's even simpler to understand potentially. What do you think about interpreting this? The real component is just the amplitude of the cosine for that frequency compared to what I was originally thinking was tapering the magnitude of the vector by it's alignment with 0 or 180 degree phase. Since those are equivalent, I'm thinking there might be a simpler interpretation, I guess that the cosine is already zero at 90 and 270 degrees so the real part already includes that tapering.

[larsoner]

I have a hard time interpreting the real or imaginary parts of these filters. In spectral analysis it means some cosine or sine part of a complex exponential (it's just a convenient way of bookkeeping). And in cross-spectrum / coherence analysis it means some phase shift between signals of interest. So plotting just the real (or imag) part of the signal only gives you part of the information. I don't know the right thing to do here but plotting just the real part (or just the imaginary part) seems to miss a lot

[alexrockhill]

Hmm those are good examples of where the real and imaginary parts of the filters are used directly but I think via the Hilbert transform, the imaginary part of the signal is just the real signal shifted in phase by a quarter cycle so I think that way the information is actually redundant and you might not be missing anything after all. In this gist (https://gist.github.com/alexrockhill/58c4316415c54fcc24f85a617862e2bf) I simulate a pure sine wave and then plot it via this method and similarly in the code below, you what I would expect to recover; a vector oscillating in the direction that was simulated. If there is a reason why visualizing it this way and this interpretation is wrong, I'm happy to abandon the idea or table until we do it some other way that is better, but this is what I was thinking would be helpful to visualize--when you get a pure oscillation, the vector changes direction forward and backward on the brain in the direction of that oscillation, I think the math on how to get there changes the interpretation but I don't quite understand why this isn't helpful or what is wrong about this visualization.

import numpy as np
import mne
import matplotlib.pyplot as plt
info = mne.create_info(['1'], 1000, 'eeg')
epochs = mne.EpochsArray(np.sin(np.linspace(0, 10, 10001) * 2 * np.pi * 20)[None, None] * 1e-6, info)
epochs_tfr = mne.time_frequency.tfr_morlet(epochs, [20], [10], return_itc=False, average=False, output='complex')
plt.plot(epochs_tfr.data.real.flatten())
plt.plot(epochs_tfr.data.imag.flatten())

[larsoner]

the imaginary part of the signal is just the real signal shifted in phase by a quarter cycle

Agreed. Same can be said of the real and imaginary parts of the Fourier basis (cos +j sin).

I think that way the information is actually redundant

No, I don't think this is how the Hilbert transform is meant to be interpreted/thought of in practice in signal processing any more than you can think of the real and imaginary coefficients of a Fourier decomposition as being redundant (the bases are actually orthogonal). And it becomes more complicated once you're talking about it in the context of a cross special density estimate, which underlies DICS, which is where I assume this would actually be useful. I think a deeper understanding of the underlying signal processing principles at play here (ideally derived from math) should guide any public implementation here rather than simulations.

[alexrockhill]

The imaginary part isn't totally redundant, but it is similar to the phase-shifted real component in actual data. The cosine of the phase scaling the magnitude is identical to the real value in the experimental data also.

import numpy as np
import mne
import matplotlib.pyplot as plt
evoked = mne.read_evokeds(mne.datasets.sample.data_path() /  'MEG' / 'sample' / 'sample_audvis-ave.fif')[0]
evoked_tfr = mne.time_frequency.tfr_morlet(evoked, [20], [2], return_itc=False, average=False, output='complex')
n = evoked.info['sfreq'] / 20 / 4
plt.plot(evoked_tfr.data[0, 0, 0].real)
plt.plot(evoked_tfr.data[0, 0, 0, int(round(n)):].imag)
plt.show()

np.allclose(evoked_tfr.data.real, np.abs(evoked_tfr.data) * np.cos(np.angle(evoked_tfr.data)))

As I understand it, the cross-spectral density is just for whitening and the complex data itself is what gets input to the minimum norm or beamformer solutions.

I don't want to post hoc justify that the real component is what is the right thing to visualize when the phase-scaled magnitude was what I was going for and was theoretically motivated but it would also be nice to visualize the direction. I know a first-component-of-the-SVD method has also been proposed which would be interesting to see implemented but it would also be nice to plot something that isn't just the direction of maximum variance if possible. Since this gets a bit complicated to understand and interpret, it would be nice if there was something with minimal processing of the data. The real component is as minimal as it gets but obviously might be too simple and not the right solution.

[alexrockhill]

Thanks for a good discussion, how about we table it for now and come back to it later?

[larsoner]

Yeah we can discuss it later