[ENH] Add option to store and return TFR taper weights

[tsbinns]

Reference issue (if any)

PR for #12851

What does this implement/fix?

Adds an option to return taper weights for complex and phase outputs of the multitaper method in tfr_array_multitaper(), and also ensures taper weights are stored in TFR objects.

Additional information

When working on this, I discovered a couple of other issues with the per-taper TFR implementations (#12851 (comment)), including the fact that the TFR object plotting methods and to_data_frame methods do not account for a taper dimension, leading to errors. Wasn't sure if people want me to also address these here or in a separate PR.

[tsbinns]

This I am somewhat unsure on. The existing implementation is to just use conc as-is, however in the MNE-Connectivity implementation that sqrt is taken: https://github.com/mne-tools/mne-connectivity/blob/97147a57eefb36a5c9680e539fdc6343a1183f20/mne_connectivity/spectral/time.py#L825

[drammock]

I am also unsure on this point. We should ask @ruuskas (who wrote the implementation in MNE-Connectivity) and @larsoner (who wrote the SciPy DPSS implementation) to weigh in.

[tsbinns]

I'm also somewhat confused about the design of the _make_dpss function:

mne-python/mne/time_frequency/tfr.py

Lines 285 to 315 in 82fc2f7

	for m in range(n_taps):
	Wm = list()
	Cm = list()
	for k, f in enumerate(freqs):
	if len(n_cycles) != 1:
	this_n_cycles = n_cycles[k]
	else:
	this_n_cycles = n_cycles[0]
	t_win = this_n_cycles / float(f)
	t = np.arange(0.0, t_win, 1.0 / sfreq)
	# Making sure wavelets are centered before tapering
	oscillation = np.exp(2.0 * 1j * np.pi * f * (t - t_win / 2.0))
	# Get dpss tapers
	tapers, conc = dpss_windows(
	t.shape[0], time_bandwidth / 2.0, n_taps, sym=False
	)
	Wk = oscillation * tapers[m]
	if zero_mean: # to make it zero mean
	real_offset = Wk.mean()
	Wk -= real_offset
	Wk /= np.sqrt(0.5) * np.linalg.norm(Wk.ravel())
	Ck = np.sqrt(conc[m])
	Wm.append(Wk)
	Cm.append(Ck)
	Ws.append(Wm)
	Cs.append(Cm)

It is looping over tapers, and then over frequencies. However, the dpss_windows function it calls internally provides the tapers and their weights for all tapers of a given frequency.

Would it not be more efficient to only loop over frequencies and take advantage of the fact that this will also return information for each taper?

[tsbinns]

I also have a question regarding testing: for the I/O tests, we're reading TFR objects that do not have a weights property (just gets assigned to None) when loaded. Do I need to create new TFR objects that actually have some weights, or is the current test sufficient?

Apart from this there are still some tests I need to expand.

[drammock]

I am also unsure on this point. We should ask @ruuskas (who wrote the implementation in MNE-Connectivity) and @larsoner (who wrote the SciPy DPSS implementation) to weigh in.

[tsbinns]

Thanks for the review @drammock! I will sort out those remaining tests, although I'm in the process of moving at the moment so it might not be for some days.

Regarding those issues I came across with TFR multitapers and converting to dataframes / plotting: would you like me to incorporate that into this PR?

[tsbinns]

Sorry for the lack of work on this, had to organise things for my PhD defence.

Everything new added here has test coverage now.

@drammock, just a couple points I would appreciate your input on:

• For the I/O tests, we're reading TFR objects that do not have a weights property (just gets assigned to None) when loaded. Do I need to create new TFR objects that actually have some weights, or is the current test sufficient?
• Regarding those issues I came across with TFR multitapers and converting to dataframes / plotting (Store n_cycles and time_bandwidth params in *TFR objects #12851 (comment)): would you like me to open a separate PR to keep individual diffs small, or just add the changes here?

Also tagging @larsoner and @ruuskas in case they can help clarify an outstanding point: #12910 (comment)

[tsbinns]

Currently working on support for a tapers dimension in ...TFRArray objects.

[drammock]

Do I need to create new TFR objects that actually have some weights, or is the current test sufficient?

Yes I think we should. most (all?) of them are created by pytest fixtures at present. I see 3 options:

1. tweak the fixtures to always return TFRs that have weights.
2. when you want to test something specific to weights, monkey-patch some weights (and a taper dim) into the object at the start of the test
3. write a new fixture (or parametrize an existing one) so that you can get TFRs with/without weights at need.

To really test thoroughly, option (2) is probably best, because then you can also patch in things that are expected to fail, and test that they do fail in the expected way.

[tsbinns]

Have removed dims being set in BaseTFR since the possibility of the optional epoch and taper dimensions makes it really difficult to disentangle here. It's much easier to handle this in the individual RawTFR, EpochsTFR, and AverageTFR classes.

[tsbinns]

Example of handling dims in the AverageTFR class where only one dimension (taper) is optional.

[tsbinns]

In terms of averaging for data with tapers, I went for the same approach we're using for Spectrum and just disallowing this.

I don't think this is an API change requiring a deprecation cycle since:

1. the docstring expects the data to not have a taper dimension, e.g. If callable, must take a NumPy array of shape (n_epochs, n_channels, n_freqs, n_times).
2. trying to call this method on an object with a taper dimension would raise an uncaught error: n_epochs, n_channels, n_freqs, n_times = self.data.shape (wouldn't be able to unpack this properly).

So explicitly preventing this method being called with a taper dimension doesn't change current behaviour, it just gives a nicer error as to why this can't be done.

[tsbinns]

It's a similar case to averaging for the time_frequency.combine_tfr() function (which also gets called by the grand_average() function).

However, unlike the EpochsTFR.average() method, this could be considered an API change since combine_tfr() should currently run with taper data. Does preventing this use case require a deprecation cycle?

On a side note, I noticed that while a public function, combine_tfr() is not listed in the API (the equivalent combine_evoked() is). Is this an oversight or an intended omission?