[ENH] Add tutorial on time-frequency source estimation with STC viewer GUI

[alexrockhill]

Related to #10721, GSoC.

So I realize there aren't really any tutorials on time-frequency source space resolution (there is DICS but that's just frequency in the tutorial) so I think this fills a very important gap in MNE functionality. Also, since it is the most complete use case (you can always average across time or trials) it covers a lot of other use cases so I'm pretty excited about adding this.

Basically, time-frequency resolved MEG data is filtered using an LCMV beamformer. This yields a list of lists of stcs where the outer layer is frequencies and the inner layer is trials. The above referenced issue is the discussion for why this is hard to keep track of and could use the introduction of a new object. There's a lot that could be done but the priorities in order are:

• Make the user interface (in pyvistaqt first)
  • Once [ENH, MRG] Add object-oriented widget abstraction #10913 is merged, use the abstraction
• Make a TFRVolSourceEstimate class to hold the list of list stcs and be returned by apply_lcmv_epochs_tfr
• Allow other functions to return corresponding TFR source objects as much as possible? This is a bit scope creep but it would be inconsistent to only have this functionality for LCMV. Making a new class for each existing source estimate class seems like a ton of overhead but potentially the more flexible source vertices+time+frequency+trials stc classes could supersede the source vertices+time source estimates and those could just be the flat case of the more general source estimates, that might be a huge refactoring though since it's so central to MNE.

Now is a great time to take five minutes to review the future direction if you have a minute. I will be out of town for the next few days so there won't be much movement.

cc @britta-wstnr @larsoner @drammock

[britta-wstnr]

I think there is some mixing of two concepts here, let me try and make this clearer:

Very generally speaking:

1. LCMV beamformers are for time-domain data (i.e., raw data, epochs, or evoked data).
2. DICS beamformers are for frequency-domain data (i.e., spectral power - in the supported case non-time-resolved power).

Some exceptions that come to mind (@alexrockhill and I have discussed some of these, so I am adding them for clarity):

1. The Hilbert beamformer or any type of bandpass-filtered signal. This is time-domain data, but (implicitly) also frequency-linked due to the filtering. But there is no power estimate (amplitude instead).
2. Our DICS code has e.g. apply_dics_evoked, which applies a DICS (i.e. frequency-resolved) beamformer to time-domain data. The use cases are not very typical I think - I would guess that people would maybe use a bandpass-filtered LCMV for this.

Now back to the proposal. I would suggest to have a DICS beamformer with this functionality. Meaning: We use make_dics to make a DICS beamformer (using a cross-spectral density matrix based on the whole epoch) and then apply this step-wise to a time-resolved power estimate (i.e., time-frequency representation, TFR).
To me, it is not so ideal to use an LCMV beamformer for that. The output would be one stc per time point we estimated power in the TFR at.

Alternatively, we can go with the Hilbert approach, where we would use bandpass-filtered time-domain data. Here the output would be one stc per frequency-band we filtered the epochs at.

Hope this will add some clarity with respect to frequency- and time-domain beamforming.

[alexrockhill]

Ok, I'm totally on the same page with one thing that I don't quite understand. I will make it more of a priority to add the DICS support sooner rather than later. The one thing I don't understand:

To me, it is not so ideal to use an LCMV beamformer for that.

Why is it appropriate to use the Hilbert transform and not any other time-frequency decomposition? You can get amplitude and not power from Morlet wavelets or you can square the Hilbert amplitude and get power right? I'm not sure what's different about the Hilbert transform compared to other time-frequency methods. From our earlier conversation, I thought you could apply a LCMV or DICS filter with the tradeoff being that the DICS is discontinuous in time and the LCMV is discontinuous in frequencies so from that I would think it would be nice to be able to use either depending on your analysis.

[britta-wstnr]

Let me be a bit more specific, hope this will make it clearer.
It is not wrong to apply a DICS filter to time-domain data or an LCMV filter to frequency-domain data. But I think it is not ideal and I am -1 on actually supporting this.

My main reason:
If you take a time-domain signal and convert it from broadband to frequency-bandpassed, you use a (frequency) filter that has different characteristics (e.g. passband) than if you take this signal and transform it into the frequency-domain (e.g. by using a Fourier transform). Ultimately, you will end up with slightly (or even really) different frequencies being resolved in the signal.
So why would you want to use e.g. an LCMV beamformer created on a bandpass-filtered signal resolving 14-17 Hz (toy example), to a power estimate that is 15 +/- 1 Hz, if you can instead apply a DICS beamformer that matches this frequency resolution perfectly?
I am not convinced this is a good idea.

Thus, for the Hilbert-case (which operates on a bandpass-filtered time-domain signal), we use an LCMV (which uses a covariance matrix based on the same bandpass-filtered signal), while for a TFR reconstruction, we should rather use a DICS beamformer that is based on a cross-spectral density matrix that uses power estimates of the same frequency resolution.

Furthermore, this way we do not have to worry about some of the not-shared features between LCMV and DICS (e.g. complex valued filters).
Plus, not supporting this might lead to less confusion in users :)

From our earlier conversation, I thought you could apply a LCMV or DICS filter with the tradeoff being that the DICS is discontinuous in time and the LCMV is discontinuous in frequencies so from that I would think it would be nice to be able to use either depending on your analysis.

Remember that the filter as such is never time-resolved (not to be mixed up with time-domain!), the time-resolution is only recovered in the application of the filter to the data (where, theoretically, it is possible to apply both LCMV and DICS filters to a time-resolved signal).

[alexrockhill]

So why would you want to use e.g. an LCMV beamformer created on a bandpass-filtered signal resolving 14-17 Hz (toy example), to a power estimate that is 15 +/- 1 Hz, if you can instead apply a DICS beamformer that matches this frequency resolution perfectly?

Because the DICS cannot be computed on a single time point and instead requires a window, so if you needed time resolution at your sampling rate but didn't mind if your frequencies were a bit smeared, I think it would be reasonable to use LCMV.

Thus, for the Hilbert-case (which operates on a bandpass-filtered time-domain signal), we use an LCMV (which uses a covariance matrix based on the same bandpass-filtered signal), while for a TFR reconstruction, we should rather use a DICS beamformer that is based on a cross-spectral density matrix that uses power estimates of the same frequency resolution.

I don't follow the rest of the argument about why a collection of bandpassed Hilbert transformed time-frequency data is any different than any other time-frequency decomposition. I understand everything you said and I think that's all reasonable but I don't understand how it applies to this point. You could take the covariance and not cross-spectral density of the TFR epochs and I'm not convinced that's any different than doing that for the Hilbert method.

Thanks for writing back, sorry for a few more questions :)

[britta-wstnr]

Because the DICS cannot be computed on a single time point and instead requires a window, so if you needed time resolution at your sampling rate but didn't mind if your frequencies were a bit smeared, I think it would be reasonable to use LCMV.

The DICS and the LCMV spatial filter are not different from each other in terms of how much data they need to go into the CSD or covariance matrix (i.e. the time interval the filter is built on).
The signal you apply the spatial filter to is of course bound in its time resolution. However, if you want the high time resolution of a Hilbert signal, why not use an LCMV beamformer for that? Why are you so focused on mixing Hilbert input with a DICS spatial filter?
If you instead use an LCMV, your frequencies in the spatial filter and the input signal you apply it to will match up perfectly.

Basically, it boils down to: why would you mix-and-match if staying within the time-domain or the frequency-domain will give you a methodologically more sound and cleaner result?

[alexrockhill]

Basically, it boils down to: why would you mix-and-match if staying within the time-domain or the frequency-domain will give you a methodologically more sound and cleaner result?

I'm still not following, you have to mix and match because you want time-frequency resolved data so you can either 1) convert to time-frequency with many bandpass filtered signals and apply Hilbert or use Morlet wavelets and then compute the covariance matrix and apply the LCMV beamformer or 2) compute the DICS beamformer on sliding time-windows using cross-spectral density, right?

To me the bandpass filter and sliding time windows are unideal because they inexactly resolve in frequency and time respectively. So wouldn't you want to select the more reasonable one for your use-case or compare results?

EDIT: After re-reading, it sounds like you associate the Hilbert transform with the time domain in some way that is different than Morlet wavelets. I would question this assumption: they're different methods but fundamentally they are trying to compute the same thing; the time-frequency amplitude or power. I've read and understand what you're saying but I don't agree that the Hilbert transform is somehow more closely associated with the time domain and therefore more cleanly associated with using an LCMV beamformer. Fundamentally, I think it's an empirical question of which works better based on that method's tradeoffs with estimating time-frequencies (i.e. choosing things like the bandwidth of the filter for Hilbert and the number of cycles for Morlet wavelets). Undoubtably, the practical experience of someone who uses these methods is very valuable and should be used to inform which method gets used or highlighted in the MNE tutorial, but I think it's important not to conflate what is mathematically inconsistent with what empirically just doesn't work as well. As I mentioned in the meeting, I think something that would really improve the methods in M/EEG in general, not just in MNE, is a really good set of decoding datasets to test methodological choices on. Maybe some of the mne-bids pipeline data that @agramfort is working on would be good for that but it's a bit outside of GSoC. Yes, if some choice of method works well for one dataset that's not super helpful in all cases, but if some methodological choice is consistently better across several datasets, I think that's much better evidence than intuition based on prior experience, again, although that is super helpful too. I don't want to dwell too much on beamforming choices but this is super relevant to the API; should the LCMV beamformer apply method take in an EpochsTFR object or Epochs and a set of frequencies and bandwidths to apply the Hilbert? I would argue that the first option is both more elegant/fits within existing MNE structures better and, until I see data or a convincing argument about why it's incoherent to use Morlet wavelets/Welch/Stockwell/whatever with LCMV, it is agnostic to the way the EpochsTFR was constructed and so is maximally abstract and generalizable. This is also on theme with trying to make a TFRVolSourceEstimate that can be returned both by the LCMV apply method, the DICS method and the apply_inverse method; although I don't want to put a lot of effort into developing all these methods, the API won't be very good if it isn't ready for these methods to be introduced in the future. That way we can highlight what is empirically shown to be the best method in the tutorial but still allow for other methods in a way that is really good for the code base because it leaves the time-frequency stuff to be done in a different, separate module of the code.

[alexrockhill]

Ok @britta-wstnr and @agramfort, if you want to try running through this tutorial, I think now would give a good indication of the direction I think we should go (and that we've agreed upon).

The main reason I bring this up is that the LCMV (of Hilbert data) and DICS apply functions returns a list of list of stcs but the LCMV is nested frequencies first and then epochs and then the stcs are shape (n_verts, n_samples) whereas the DICS is nested sliding time windows first then epochs and then the stcs are shape (n_verts, n_freqs). I think it's pretty difficult to rearrange the dimensions in lists of list so it really seems to be crying out for a class abstraction that will arrange them correctly for you. That's the next step that I'll try and get done today which is the main thing holding up making an actual GUI.

As for whether it's okay to use Morlet wavelets instead of Hilbert and why you can't use power as input to a minimum norm estimation, I think we should discuss at the next meeting and I'm looking forward to learning about how these work :)

It's not crucial that everything be resolved with the source estimation before moving forward (it really needs to be optimized for parallelization also), I just wanted to make sure the class worked so that it could be used as input to the GUI without having to be changed.

I have no idea what the data looks like and I haven't invented any new methods, just used the old ones, so I assume that it will be reasonable but a time-frequency viewer is sorely needed with this high of a dimension of data; it's just very hard to look at. Will have a reasonable screenshot soon!

EDIT: Also very likely the final minimum norm source estimation method is wrong since it's done with wavelets only in mne.minimum_norm.source_induced_power but it would be good to know why that is so. I'll switch it over to that code actually now that I think about it.

EDIT2: Actually I tried source_induced_power and waited five minutes but it didn't finish computing so I'll wait to discuss the best options.

[alexrockhill]

Notes per method:

1a) The LCMV on the Hilbert transformed data works great but is quite slow because the bandpass filter has to be computed on the raw data before epoching.
1b) The LCMV on Morlet wavelets seems to be efficient and has the fewest drawbacks computationally but we'll see about the results
2) The DICS is quite slow as well and there are issues with the sliding window size for lower frequencies since it must be larger. I had to go up to 90 Hz in the tests with the testing sampling frequency of 300 Hz because I kept getting warnings that the filter was longer than the data. This seems like a big drawback for this method in making a single TFR source estimate with the same sliding windows. If you allowed for bigger and smaller windows the code would be very complicated as it would no longer be matrix based... this is a tough one.
3) I imagine just applying the minimum norm to time-frequency epochs data may be suboptimal and incorrect but it computes in a reasonable time everything works well with window sizes since the time-frequencies are computed beforehand.

[agramfort]

@alexrockhill it would help to go over a rendered doc page / tutorial to see the user perspective on this. Is there any? Sorry if I let you babysit me here ...

[alexrockhill]

@alexrockhill it would help to go over a rendered doc page / tutorial to see the user perspective on this. Is there any? Sorry if I let you babysit me here ...

There will be soon since #10940 merged, it was too much memory, causing a segfault in circle but it runs locally. Let me see about now, otherwise I'll comment out 2/3 methods (LCMV, DICS, MNE).

[britta-wstnr]

Okay, I can only reiterate my points from before:

1. Hilbert outputs the same time-resolution as you put in, thus I sort this into a time-domain signal. The Morlet wavelet approach mostly does not do this - you traditionally have not the same time steps as in your input signal (but fewer).
2. And more importantly: The Hilbert uses a time-domain filter. The covariance matrix of the LCMV beamformer can utilize a time domain filter. The Morlet wavelet uses a wavelet. The CSD matrix of the DICS beamformer uses a wavelet. Thus, staying "within" domain makes the beamformer filter and the signal you apply the filter to match way better.

I am still opposed to applying an LCMV beamformer to a Morlet wavelet resolved TFR. Of course, you could show this actually works way better than the DICS beamformer on a Morlet TFR - but that would be a whole project (you cannot just show it on one data set and be done but need carefully simulated data etc).

Until then: let's not confuse our users.

[alexrockhill]

1. Hilbert outputs the same time-resolution as you put in, thus I sort this into a time-domain signal. The Morlet wavelet approach mostly does not do this - you traditionally have not the same time steps as in your input signal (but fewer).

That's not my understanding, here is the API, https://mne.tools/stable/generated/mne.time_frequency.tfr_array_morlet.html, you can decimate but by default there are the same number of time steps. And you can decimate the Hilbert in the same way.

2. And more importantly: The Hilbert uses a time-domain filter. The covariance matrix of the LCMV beamformer can utilize a time domain filter. The Morlet wavelet uses a wavelet. The CSD matrix of the DICS beamformer uses a wavelet. Thus, staying "within" domain makes the beamformer filter and the signal you apply the filter to match way better.

I agree that the CSD matrix uses a wavelet and thus matching CSD methods with your filter and and time-series object that you apply it to seems like the only logically consistent approach but I don't know what this has to do with LCMV, that seems like a misleading analogy. The Hilbert transform itself doesn't use a time-domain filter, you first filter your data before taking the analytic representation of the data which is very similar to a Fourier transform. Both Hilbert and Morlet methods generally window data in frequency (via the filter before Hilbert) and time-frequency respectively and do a Fourier transform or Fourier-transform-like operation to compute the amplitude or power so I don't see the big difference.

Happy to talk more in person, I do think by iterating we've clarified where each of our thought processes are coming from but it seems like I haven't done a good job convincing you so maybe we should see if other people weighing in helps us get on the same page.

EDIT: Maybe @wmvanvliet might have a minute to chime in here since you've done a lot of the DICS development. The question is: Is it okay to use Morlet wavelets instead of narrow bandpass filtered epoch data with a Hilbert transform as input to an LCMV beamformer? @britta-wstnr 's and my points of view is summarized above, if you have a few minutes.

[alexrockhill]

Here's the artifact @agramfort https://output.circle-artifacts.com/output/job/10323476-dfd6-4c9a-86c0-58228c1aa59f/artifacts/0/dev/auto_tutorials/inverse/55_tfr_stc.html?highlight=55_tfr, sorry I know it's late your time. It doesn't have many visuals because that has to be done when the GUI is built but basically:

1. You narrowband bandpass filter the data and apply the Hilbert transform to make a numpy array and then use that to construct an EpochsTFR object from scratch, then you compute LCMV filters and apply them which are now returned as an stc object with dimensions (n_epochs, n_verts, n_orient, n_freqs, n_times).

2. You compute cross-spectral density on sliding time windows to get the same stc object with DICS but this time the times are the center of the windows so there are fewer of them but the frequency resolution is really good.

3. You do the same thing as LCMV but treat the time-frequency data just like time series data and apply standard MNE looping over frequencies.

[alexrockhill]

Fixes #10721.

I only mentioned passingly in a comment to a review from @britta-wstnr that I would like to keep the epochs dimension for my own analysis at single-trial resolution if that's not too much of a headache for everyone else.

Also, I didn't inherit TFRSourceEstimate from _BaseSourceEstimate yet because I thought it might be nice to use the newer info structure and inherit the mixins such as those being refactored here #10945. I'm happy to inherit if that's too iconoclastic but I thought it might be a good time to start that kind of refactoring and build everything around mne.viz.Brain and the GUI abstraction #10913. I won't be able to refactor all the other MNE SourceEstimates this summer but I think it's not too redundant to maintain during a transition period.

[alexrockhill]

@agramfort and @britta-wstnr , I get the following stack trace when I try to use complex phase-amplitude and not power:

ValueError                                Traceback (most recent call last)
<ipython-input-2-a143c40a3035> in <cell line: 11>()
     14     data_rank.update({ch_type: ch_rank for ch_type, ch_rank in
     15                       noise_rank.items() if ch_rank < data_rank[ch_type]})
---> 16     filters.append(mne.beamformer.make_lcmv(
     17         epochs_tfr.info, forward, data_cov, reg=0.05,
     18         noise_cov=noise_cov, pick_ori='vector',

<decorator-gen-497> in make_lcmv(info, forward, data_cov, reg, noise_cov, label, pick_ori, rank, weight_norm, reduce_rank, depth, inversion, verbose)

~/projects/mne-python/mne/beamformer/_lcmv.py in make_lcmv(***failed resolving arguments***)
    174     # compute spatial filter
    175     n_orient = 3 if is_free_ori else 1
--> 176     W, max_power_ori = _compute_beamformer(
    177         G, Cm, reg, n_orient, weight_norm, pick_ori, reduce_rank, rank_int,
    178         inversion=inversion, nn=nn, orient_std=orient_std,

~/projects/mne-python/mne/beamformer/_compute_beamformer.py in _compute_beamformer(***failed resolving arguments***)
    310     # Here W is W_ug, i.e.:
    311     # G.T @ Cm_inv / (G.T @ Cm_inv @ G)
--> 312     bf_denom_inv = _sym_inv_sm(bf_denom, reduce_rank, inversion, sk)
    313     assert bf_denom_inv.shape == (n_sources, n_orient, n_orient)
    314     W = np.matmul(bf_denom_inv, bf_numer)

~/projects/mne-python/mne/beamformer/_compute_beamformer.py in _sym_inv_sm(x, reduce_rank, inversion, sk)
    120         assert x.shape[1:] == (3, 3)
    121         if inversion == 'matrix':
--> 122             x_inv = _sym_mat_pow(x, -1, reduce_rank=reduce_rank)
    123             # Reapply source covariance after inversion
    124             x_inv *= sk[:, :, np.newaxis]

~/projects/mne-python/mne/utils/linalg.py in _sym_mat_pow(A, power, rcond, reduce_rank, return_s)
    170     limit = s[..., -1:] * rcond
    171     if not (s >= -limit).all():  # allow some tiny small negative ones
--> 172         raise ValueError('Matrix is not positive semi-definite')
    173     s[s <= limit] = np.inf if power < 0 else 0
    174     if reduce_rank:

ValueError: Matrix is not positive semi-definite

Any thoughts on how to fix this?

EDIT: I'm pushing a hacked version for now that uses single dipole inversion and suppressed an error so that you can see the changes but it definitely needs to be fixed.

The minimum norm source estimation worked with no errors also by the way.

[alexrockhill]

Ok, I couldn't figure out why the eigenvalues were not all positive but it has to do with the data covariance matrix which is positive definite for the CSD and is complex as well in that case so complex numbers are supported, it's just an issue with the conditioning of the matrix. The CSD matrix had a much larger diagonal than off diagonal entries for the real part compared to the covariance of the Hilbert which had larger off the diagonal than on the diagonal. I'm not exactly sure what the issue is with the conditioning though, I'm a bit stuck.

As I'm writing this, I'm thinking that real part of the covariance of the Hilbert should be positive since it's amplitude and I didn't check that this was true, I'll look into that next.

[agramfort]

hum hard to say without debugging. I think the code uses eigh which should also work with a complex valued matrix. Did you see how negative the eigenvalues are? Can you check at least that the input matrix is hermitian is A == A.H ?

[wmvanvliet]

Maybe @wmvanvliet might have a minute to chime in here since you've done a lot of the DICS development. The question is: Is it okay to use Morlet wavelets instead of narrow bandpass filtered epoch data with a Hilbert transform as input to an LCMV beamformer?

DICS == LCMV applied to a TFR (wavelets, multitaper, hilbert, whatever). A "tfr covariance" == CSD. So if you do LCMV using a "covariance matrix" derived from some sort of time-frequency decomposition, you are doing DICS. Literally, DICS and LCMV use the same code path (_compute_beamformer) with the only difference being s/cov/csd. So, it makes no sense to have an example that does this "manually" without using the DICS code.

What you are doing here is event-related DICS (Laaksonen, Kujala and Salmelin 2007) which we never got around to implementing. It looks like you are doing this with the apply_dics_evoked_tfr function, which is greatly appreciated.

[alexrockhill]

That was my understanding but now that you say that maybe it doesn't make sense to have the sliding windows but rather to only use the compute_covariance function on the epochs_tfr.

[alexrockhill]

hum hard to say without debugging. I think the code uses eigh which should also work with a complex valued matrix. Did you see how negative the eigenvalues are? Can you check at least that the input matrix is hermitian is A == A.H ?

They are large and negative, the largest was -13000 or so. I'll have to check on Hermitian.

[wmvanvliet]

now that you say that maybe it doesn't make sense to have the sliding windows

I think it is best to follow published methods for the main MNE-Python package. So I recommend following whatever approach they used in the paper. A sliding window will always have to be applied somewhere. When computing the TFR, when computing the CSD, or when applying the beamformer... I guess in the end it doesn't matter where we apply it, or even if we apply multiple sliding windows.

[wmvanvliet]

the compute_covariance function on the epochs_tfr.

Could you please modify the compute_csd function to be able to operate on TFRs instead? Computing "covariance" on a TFR object means computing a CSD.

[britta-wstnr]

I agree that the CSD matrix uses a wavelet and thus matching CSD methods with your filter and and time-series object that you apply it to seems like the only logically consistent approach but I don't know what this has to do with LCMV, that seems like a misleading analogy.

The LCMV beamformer uses a covariance matrix based on time-domain data. The Hilbert beamformer, as it was defined by my former supervisor @sarangnemo, computes a covariance matrix based on filtered time-domain data and then uses this in the LCMV beamformer.
If you would now use a DICS for the same operation instead, you would need a CSD matrix as input, e.g. obtained via Morlet wavelets (as implemented in MNE), and the exact frequencies resolved would not match anymore. Especially so, since a Hilbert beamformer normally uses a broader frequency band (to achieve higher time resolution/less smearing in the time domain per the time-frequency tradeoff), which is hard to achieve with wavelets or Fourier (you'd need multitaper for that). That mismatch in frequencies resolved is what I am suspecting to be a really bad idea in theory and practice, what has not been done before in practice, and what I see as misleading for the user (as many will not be aware of the mismatch in frequency resolution).

I would vote for not making the focus to invent a new beamformer, but keep with what is used in the literature already - and for your GSOC just use synthetic input as suggested by @agramfort. The details of the computation of the beamformer should not influence the layout of the source space object.
As mentioned before, in cases where it would make more sense to compute the spatial filter iteratively (Hilbert case, because you also have to do the filtering etc, which might make more sense to do on the fly in a loop instead of first creating all those epoch objects in different frequencies and with and without Hilbert applied), the output could still be transformed into a multi-dimensional source space object after the fact.

[britta-wstnr]

DICS == LCMV applied to a TFR (wavelets, multitaper, hilbert, whatever). A "tfr covariance" == CSD. So if you do LCMV using a "covariance matrix" derived from some sort of time-frequency decomposition, you are doing DICS. Literally, DICS and LCMV use the same code path (_compute_beamformer) with the only difference being s/cov/csd. So, it makes no sense to have an example that does this "manually" without using the DICS code.

What I would do is computed on CSD matrix across power for the whole time period - and then apply this time-step by time-step --- as you would with an LCMV beamformer that you apply to an evoked.
For everything else we get very close to the 5-D beamformer by @sarangnemo, which we have just deprecated (as it is just not the best way of doing things anymore: your spatial filter suffers under not having many samples for your covariance/CSD matrix if you "slide" the make_filter operation across time.

[britta-wstnr]

I will try and check out your changes and the behavior of the GUI.
There seem to be a lot of changes as per your latest comments, and I am not sure what they address - are they directly related to reviews or are they adding new functionality?
It's hard for me to see what changed because the "View changes since last review" is broken for me here (probably due to rebase or force-push?) and there is a lot of commits and lines of code to check out.
Thus, this will take me some time to get through - could you share a brief list of what functionalities you changed/added (if any) beyond the lists/generator and cosmetic fixes we had asked for?

[agramfort]

@alexrockhill here is a batch comments. I can certainly do another pass next week end. Can you have rebasing your PRs? It forces me to do reset --hard and we squash merge PRs anyway.

Regarding UI

It seems the min is not used over the MRI slice. Any voxel that is less that min should be transparent for me so I see the MRI. Here it's always black so I never see the MRI.

🙏

[agramfort]

one reason I would rather use csd_ers rather than ers_csd is that when I type I prefer to autocomplete based on the type of object. Basically all variables that are of type csd start with csd_ so csd_ gives me what I want.

[alexrockhill]

Sounds good, no problem, again, I was just changing it to be consistent, there was csd_baseline and active_cov etc. it would just be nice if the one that was first was consistent. Again, not trying to change all of MNE but just to be consistent in this tutorial would be nice. I will change those to the type first and the descriptor second like you're asking.

[agramfort]

where is this 1e16 coming from?

[alexrockhill]

If you cast to integers without rescaling, you mostly get zeros so this is just so that you can see your data but use integers for speed like is often done for MRIs.

[alexrockhill]

Ok, I will stop rebasing for merge conflicts. @britta-wstnr, the changes were from the request that @drammock and @agramfort had that the function take in the list of lists of stcs directly. If that is done without modification, it's unusable on my computer at least, it's far too slow. So I just changed to put the casting to integers in the function. That's where the 1e16 etc comes from @agramfort also.

Making the suggested changes now.

Also, I empathize with the added overhead for reviewing a PR but I wouldn't have touched the comments if they weren't inconsistent within that tutorial (some inline comments say Compute the TFR here. and others say compute the TFR here) this to me is added overhead for both the user trying to figure out why a tutorial seems to be written by three different people and the contributor trying to match style. I'm not trying to go through all of MNE and say that all inline comments should be informal (i.e. no capitalization) but to me, it just seems common sense to fix things when you come across them. I will revert since I'm getting so much push back on this but this seems like a no-brainer to me and I really don't understand the perspective of why not to fix this (other than of course if you take it to the extreme of fixing everything in MNE which again is not what I was trying to do).

[alexrockhill]

I tried with the typecasting to get everything faster and it was working, it just adds a lot of complexity so I thought better of it for now, maybe will revisit later. This version should pass all the tests and work well so long as you have a decently powered computer.

[alexrockhill]

I took a minute and fixed the integer casting, it's not really doable for ITC, that was a hangup, that will have to be cast to floating point numbers for the whole data so that will be a limitation for larger datasets. This version was very fast on my machine so I think it would work if there were more epochs and/or more frequencies which I think is pretty important for this scaling well.

[larsoner]

@agramfort ready for another review?

Failures are unrelated I think

[agramfort]

thx @alexrockhill for the updates !

here are a few comments on the code.

now on the GUI I don't know if it's only me but the time and frequency sliders are quickly messed up in their locations when I click on the 2d TFR image at the bottom. Are the sliders and the cross in the image allows consistent at your end?

also do you know what takes the most time when I click on TFR image it can take 3s to update. Is it for you the same? what slows the GUI? data access, copies? matplotlib?

🙏

[agramfort]

ok fair enough

[agramfort]

in the header docstring above you say

"""View a volume time and/or frequency source time course estimate.

I therefore wondering if it should not be a and only.
Do you have a example of usage without frequency?

[alexrockhill]

Yes, it works just fine without frequency

[alexrockhill]

There is an example in the tests, do you want to add it to a tutorial or example?

[agramfort]

would there be a way to avoid having some method logic in the GUI? basically we could have other stuff than ITC, or average_power (think real part of the coherence etc.) and I am a bit worried to see this implemented in a GUI code

just tell me what you think

[alexrockhill]

The reason to implement it here is that it's really slow using float and it's relatively easy and much faster to use integers. There's not a cleanly abstracted function for ITC in tfr.py anyway so that would need to be refactored. We could make _itc and _itc_int functions in tfr.py if you want, I just was trying to keep the number of files touched low.

[alexrockhill]

It doesn't look like the most straightforward abstraction/refactor from here since it's looped

# Inter-trial phase locking is apparently computed per taper...
if 'itc' in output:
    plf = np.zeros((n_freqs, n_times), dtype=np.complex128)

# Loop across epochs
for epoch_idx, tfr in enumerate(coefs):
    # Transform complex values
    if output in ['power', 'avg_power']:
        tfr = (tfr * tfr.conj()).real  # power
    elif output == 'phase':
        tfr = np.angle(tfr)
    elif output == 'avg_power_itc':
        tfr_abs = np.abs(tfr)
        plf += tfr / tfr_abs  # phase
        tfr = tfr_abs ** 2  # power
    elif output == 'itc':
        plf += tfr / np.abs(tfr)  # phase
        continue  # not need to stack anything else than plf

    # Stack or add
    if ('avg_' in output) or ('itc' in output):
        tfrs += tfr
    elif output in ['complex', 'phase'] and method == 'multitaper':
        tfrs[taper_idx, epoch_idx] += tfr
    else:
        tfrs[epoch_idx] += tfr

[agramfort]

same remark we have some shared functions to compute baseline in different ways. I would try to rely on these has much as possible

[alexrockhill]

I looked into these, they require a time range and I just hadn't implemented it yet. I should just add a tmin and tmax for the baseline

[alexrockhill]

now on the GUI I don't know if it's only me but the time and frequency sliders are quickly messed up in their locations when I click on the 2d TFR image at the bottom. Are the sliders and the cross in the image allows consistent at your end?

Hmm, these seem like they're in sync to me. Let me know if this is still an issue.

also do you know what takes the most time when I click on TFR image it can take 3s to update. Is it for you the same? what slows the GUI? data access, copies? matplotlib?

It was that the time was accessing data and updating the images and then updating frequency did the same. I fixed it so that it only updated once.

[alexrockhill]

Copying a minimally reproducible script to run the viewer:

import numpy as np
import mne
from mne.cov import compute_covariance
from mne.datasets import somato
from mne.time_frequency import csd_tfr
from mne.beamformer import (make_dics, apply_dics_csd, make_lcmv,
                            apply_lcmv_cov)
from mne.minimum_norm import (make_inverse_operator, apply_inverse_cov)

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()

# the DICS beamformer currently only supports a single sensor type,
# we'll use the gradiometers in this example
picks = mne.pick_types(raw.info, meg='grad', 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)
baseline_cov = compute_covariance(epochs, tmin=baseline_win[0],
                                  tmax=baseline_win[1], method='shrunk',
                                  rank=rank, verbose=True)
active_cov = compute_covariance(epochs, tmin=active_win[0], tmax=active_win[1],
                                method='shrunk', rank=rank, verbose=True)

# weighted averaging is already in the addition of covariance objects
common_cov = baseline_cov + active_cov

# 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.view_vol_stc(stcs, subject=subject, subjects_dir=subjects_dir,
                              src=vol_src, inst=epochs_tfr)
viewer.go_to_max()  # 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)

Looks good to me

[larsoner]

Merge reverted, feel free to respond to above comments @agramfort and then presumably @alexrockhill can open a second PR with the same changeset for any additional commits