[ENH, MRG] Add EpochsSpectrumArray and SpectrumArray classes #11803
Conversation
|
Looks like the failures are unrelated, good to go by me |
mne/time_frequency/spectrum.py
Outdated
| inst_type_str="Raw", | ||
| data_type="Average Power Spectrum", |
There was a problem hiding this comment.
If the data comes from an array we don't know if it was computed on Raw or averaged (Evoked) data. Also we don't know if it's power/amplitude, real/complex, etc.
| inst_type_str="Raw", | |
| data_type="Average Power Spectrum", | |
| inst_type_str="Array", | |
| data_type="Unknown", |
This will probably break some tests (assuming we're testing thoroughly) so hopefully those test failures can guide what needs to change elsewhere
There was a problem hiding this comment.
So real/complex can be inferred from the data type. We ask specifically for power in the constructor so maybe we should require that and let the amplitude kwarg in Spectrum.plot handle that? That'd be my 2 cents but not highly opinionated, whatever is simplest
|
Okay, I think the outstanding things here are:
Otherwise looks good, thanks for the review @drammock. Oh one more thing, I set |
mne/time_frequency/spectrum.py
Outdated
| %(drop_log)s | ||
| %(events_epochs)s | ||
| %(event_id)s | ||
| %(metadata_epochs)s |
There was a problem hiding this comment.
did metadata get missed here or did you intentionally keep it? To me that one is the clearest case of "can / should be set afterward, not in constructor"
There was a problem hiding this comment.
I left it in, I could see a use case for it here but that sounds reasonable to modularize if we prefer they set metadata later, I'll take it out
|
Okay to merge? The |
|
Ok for me
|
Your checklist above has only 1 of 3 items checked off (plus there's a fourth item that doesn't have a checkbox --- the question about |
|
Not a huge rush but I'm sharing a script that uses it with someone else so it's a bit of a pain to describe how to install a feature branch so it would be nice not to take too long. I was going to let whoever merged the PR check those boxes to make sure they're okay with them, it's fine by me how it is right now (after I implement the latest review). |
|
I just went ahead and added the fixtures so everyone wins 😄 Oh man that's a creepy smile, maybe 🙂 that's better |
|
Looks like the failure is unrelated. |
|
sorry, I got captured by a grant deadline, it was on my list to come back to this yesterday but I ran out of time. Will prioritize for today. |
|
No worries, next week is great too, just didn't want it to get forgotten. |
|
we chatted about this PR in the dev meeting today, to gauge opinions on whether supporting complex data as input to
All this has made me realize something that I previously overlooked / didn't discuss in the meeting, which is that when we originally added complex support, we were outputting numpy arrays so the ramifications didn't need to be traced through the rest of our API to make sure the complexness of the data didn't break things. When I created the Spectrum classes, support for complex unaggregated multitaper output was included for legacy reasons, but it probably would have been better to not support it in the classes (and to keep open a separate array-yielding code path for users who want the unaggregated multitaper output). Which leads me to ask @mmagnuski (who originally requested the complex unaggregated multitaper output) a few questions:
|
|
Those are pretty reasonable concerns: you have to pass the method if it's complex and I added checks that the dimensions matched the size of the frequencies passed so I think it the checks for the wrong data make it very safe (i.e. even if frequencies and |
|
@drammock |
|
what's the status here? I see some remaining comments by @drammock. It's a matter of agreeing if we support complex or not? |
Yes, that's my understanding we should agree whether to support complex psds. One last point: I think from a maintenance perspective having complex spectrumarrays was actually very helpful because in testing complete cases, there were several bugs in how plotting was handled. |
mne/time_frequency/spectrum.py
Outdated
| # handle unaggregated multitaper | ||
| if hasattr(self, "_mt_weights"): | ||
| logger.info("Aggregating multitaper estimates before plotting...") | ||
| psds = _psd_from_mt(psds, weights=self._mt_weights) | ||
| # handle unaggregated Welch | ||
| elif "segment" in self._dims: | ||
| logger.info("Aggregating Welch estimates (median) before plotting...") | ||
| seg_axis = self._dims.index("segment") | ||
| psds = np.nanmedian(psds, axis=seg_axis) | ||
| if np.iscomplexobj(psds): # convert to power for plotting | ||
| psds = (psds * psds.conj()).real | ||
| if "epoch" in self._dims: | ||
| psds = np.mean(psds, axis=self._dims.index("epoch")) |
|
I'll restate her the pros and cons from my comment above:
@mmagnuski's comment suggests to me that we should probably deprecate complex output support in the My vote here is:
|
|
That's fine with me, I'll remove complex support, you spent a lot of time and effort making spectrum classes @drammock so I'm happy to differ to your judgment and I agree there are not a lot of apparent use cases that come to mind--phase doesn't make any sense, it's not time-frequency, thanks for the direct communication, that helps move forward. |
|
Ok I removed complex support and this is ready for review. I changed the |
There was a problem hiding this comment.
here's a self-review as commentary on the changes I've just pushed.
@larsoner @agramfort or @mmagnuski maybe best for one of you to push the green button on this one, it's got a fair amount of code from both me and @alexrockhill
| """Get raw with power spectral density computed from mne.io.tests.data.""" | ||
| return [raw.compute_psd(method=method) for method in ("welch", "multitaper")] | ||
| return raw.compute_psd() |
There was a problem hiding this comment.
I don't think we have a good reason to test both Welch and multitaper in the fixture; the differences aren't relevant to most tests (and should be covered adequately by the tests of the welch/mutitaper array methods used under the hood)
| @@ -298,9 +298,9 @@ def raw_ctf(): | |||
|
|
|||
|
|
|||
| @pytest.fixture(scope="function") | |||
| def raw_psds(raw): | |||
| def raw_spectrum(raw): | |||
mne/time_frequency/__init__.py
Outdated
| from .ar import fit_iir_model_raw | ||
| from .multitaper import dpss_windows, psd_array_multitaper, tfr_array_multitaper | ||
| from .spectrum import ( | ||
| EpochsSpectrum, | ||
| EpochsSpectrumArray, | ||
| Spectrum, | ||
| SpectrumArray, | ||
| read_spectrum, | ||
| ) | ||
| from ._stft import stft, istft, stftfreq | ||
| from ._stockwell import tfr_stockwell, tfr_array_stockwell |
There was a problem hiding this comment.
when you merged in main, this wasn't properly integrated into the new lazy loading scheme
mne/time_frequency/spectrum.py
Outdated
| @@ -429,7 +429,7 @@ def _check_values(self): | |||
| assert len(self._dims) == self._data.ndim, (self._dims, self._data.ndim) | |||
| assert self._data.shape == self._shape | |||
| # negative values OK if the spectrum is really fourier coefficients | |||
| if np.iscomplexobj(self._data): | |||
| if "taper" in self._dims: | |||
There was a problem hiding this comment.
reverting to reduce churn, esp. given that we now plan to deprecate support for complex data anyway.
mne/time_frequency/spectrum.py
Outdated
| # handle unaggregated multitaper | ||
| if hasattr(self, "_mt_weights"): | ||
| logger.info("Aggregating multitaper estimates before plotting...") | ||
| psds = _psd_from_mt(psds, weights=self._mt_weights) | ||
| # handle unaggregated Welch | ||
| elif "segment" in self._dims: | ||
| logger.info("Aggregating Welch estimates (median) before plotting...") | ||
| seg_axis = self._dims.index("segment") | ||
| psds = np.nanmedian(psds, axis=seg_axis) | ||
| if np.iscomplexobj(psds): # convert to power for plotting | ||
| psds = (psds * psds.conj()).real |
mne/time_frequency/spectrum.py
Outdated
| # handle unaggregated multitaper | ||
| if hasattr(self, "_mt_weights"): | ||
| logger.info("Aggregating multitaper estimates before plotting...") | ||
| psds = _psd_from_mt(psds, weights=self._mt_weights) | ||
| # handle unaggregated Welch | ||
| elif "segment" in self._dims: | ||
| logger.info("Aggregating Welch estimates (median) before plotting...") | ||
| seg_axis = self._dims.index("segment") | ||
| psds = np.nanmedian(psds, axis=seg_axis) | ||
| if np.iscomplexobj(psds): # convert to power for plotting | ||
| psds = (psds * psds.conj()).real |
| @@ -1226,11 +1197,21 @@ def __getitem__(self, item): | |||
| return BaseRaw._getitem(self, item, return_times=False) | |||
|
|
|||
|
|
|||
| def _check_data_shape(param, dim, data, expected): | |||
| if data.shape[dim] != expected: | |||
| def _check_data_shape(data, freqs, info, ndim): | |||
There was a problem hiding this comment.
I reworked this to check everything at once rather than needing multiple calls to the check function. Seemed simpler this way once we knew we only had to handle 2D and 3D cases.
| @@ -141,18 +142,21 @@ def test_spectrum_io(inst, tmp_path, request, evoked): | |||
| assert orig == loaded | |||
|
|
|||
|
|
|||
| def test_spectrum_copy(raw): | |||
| def test_spectrum_copy(raw_spectrum): | |||
There was a problem hiding this comment.
Although I'd have preferred adding the fixture in its own PR, I decided it was expedient to just go with it, so I've used it in the other spectrum tests too.
|
All green, this one is ready for review/merge |
There was a problem hiding this comment.
Pretty short diff in the end, thanks @alexrockhill @drammock !
…ls#11803) Co-authored-by: Daniel McCloy <[email protected]>
Makes API a lot nicer and consistent with MNE related to discussion here: https://mne.discourse.group/t/instantiating-epochsspectrum-and-spectrum-class-objects/7127.
Without the docstrings of the array, you actually have to look in
__set_state__to figure out what is required in thestateargument which is pretty confusing for me. The use case of ragged/unequal epochs is something I ran into too.Okay with you @drammock ?