Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[MAINT, MRG] Refactor TimeMixin to include more time operations, ignore filter warning for decim for EpochsTFR #10945

Merged
merged 13 commits into from
Jul 27, 2022

Conversation

alexrockhill
Copy link
Contributor

@alexrockhill alexrockhill commented Jul 21, 2022

Fix leftovers from #10940.

So I thought, might as well loop in Evoked with this to be even more DRY but in doing so it turns out that they actually decimate in different ways. It comes down to the order of the following lines:

Evoked method:

self._decim *= decim
start_idx = int(round(-self._raw_times[0] * (self.info['sfreq'] *
                                                     self._decim)))

Epochs method:

start_idx = int(round(-self._raw_times[0] * (self.info['sfreq'] *
                                                     self._decim)))
self._decim *= decim

What this leads to is:

Original times:

array([-1. , -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,  0. ,
        0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9])

Evoked method:

array([-1. , -0.4,  0.2,  0.8])

Epochs method:

array([-0.6,  0.0,  0.6])

cc @agramfort @drammock @larsoner

Seems fairly consequential...

@alexrockhill alexrockhill changed the title [MAINT, WIP] Ignore filter warning for decim for EpochsTFR [MAINT, MRG] Ignore filter warning for decim for EpochsTFR Jul 21, 2022
@alexrockhill
Copy link
Contributor Author

Also there's a bit of an issue that the mixins docstring can only say one return type when self is returned which leads to issues like here https://mne.tools/stable/generated/mne.Evoked.html#mne.Evoked.shift_time where evoked says it returns an Epochs object.

@alexrockhill
Copy link
Contributor Author

After a bit of digging, it seems like it's because the Epochs must have 0 because that's the sample that they are defined based on, which seems reasonable, it just would make sense that the Evoked would do the same.

@alexrockhill
Copy link
Contributor Author

Ok this should fix all the things I found that were inconsistent.

@alexrockhill
Copy link
Contributor Author

That looks pretty sensible to me for the multiple return type https://output.circle-artifacts.com/output/job/76871b03-098a-4353-9ac0-402f230ccd3b/artifacts/0/dev/generated/mne.Evoked.html?highlight=evoked#mne.Evoked.decimate, definitely an improvement over main.

@alexrockhill
Copy link
Contributor Author

I realized I wasn't helping too much by putting another mixin that was redundant with TimeMixin so I just refactored everything within TimeMixin which really trimmed down a lot of redundant code. This got a bit more than I thought but I think this is a really nice refactoring.

@alexrockhill alexrockhill changed the title [MAINT, MRG] Ignore filter warning for decim for EpochsTFR [MAINT, MRG] Refactor TimeMixin to include more time operations, ignore filter warning for decim for EpochsTFR Jul 22, 2022
@alexrockhill
Copy link
Contributor Author

alexrockhill commented Jul 22, 2022

Ah shoot, times just has to be set with the _set_times setter one more place and the text of a warning changed, I'll fix it tomorrow morning.

@drammock
Copy link
Member

That looks pretty sensible to me for the multiple return type https://output.circle-artifacts.com/output/job/76871b03-098a-4353-9ac0-402f230ccd3b/artifacts/0/dev/generated/mne.Evoked.html?highlight=evoked#mne.Evoked.decimate, definitely an improvement over main.

I'm not sure about the glossary entry for "time-series object". Try git grep "inst :" and you'll see lots of examples of how we handle this situation elsewhere. We're not perfectly consistent in how we do it, but let's not make that worse by adding a new variant.

@drammock
Copy link
Member

cross-ref to #10947. Can you touch those two examples/tutorials in this PR, so they will render here and we can be sure they work again?

Copy link
Member

@drammock drammock left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this looks like a nice refactoring. LMK when it's ready for another look.

doc/changes/latest.inc Outdated Show resolved Hide resolved
@@ -81,7 +81,7 @@ def test_decim():
info = create_info(n_channels, sfreq, 'eeg')
with info._unlock():
info['lowpass'] = sfreq_new / float(decim)
evoked = EvokedArray(data, info, tmin=-1)
evoked = EvokedArray(data, info)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why change this? I get that alignment will change now, so the exact times will change on the assert lines below, but an EvokedArray with negative tmin should still work.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It doesn't work with the immediately following tests with a simpler [::decim] because that doesn't account for the API change

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't follow. The API change is that decim now shifts which samples it keeps in order to ensure that it keeps the sample at t=0. That shouldn't mean that there are no negative times.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you have negative times then checking that the decimation worked properly is super complicated and done in the tests for epochs so better not to do that here.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah, OK, I see now. The lines 88 and 89 were hidden from me in the GitHub PR review UI. If those lines are the issue, then it ought to be possible to revert the change to this line, and edit the starting sample (e.g., the 1 in 1::decim) in lines 88-89 right?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is that simple. lines 84-90:

    evoked = EvokedArray(data, info, tmin=-1)  # back the way it started
    evoked_dec = evoked.copy().decimate(decim)
    evoked_dec_2 = evoked.copy().decimate(decim, offset=1)
    evoked_dec_3 = evoked.decimate(dec_1).decimate(dec_2)
    assert_array_equal(evoked_dec.data, data[:, 4::decim])    # just change the 
    assert_array_equal(evoked_dec_2.data, data[:, 5::decim])  # starting samp
    assert_array_equal(evoked_dec.data, evoked_dec_3.data)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also, please stop marking conversations as "resolved" just because you disagree. Sometimes it's obvious (like "please make this an f-string") and it's fine, but in more complex cases, whether something is "resolved" is generally best determined by the person who asked the question.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

-1 from me, the numbers 4 and 5 are picked out of nowhere, it's properly tested here:

raw, events, picks = _get_data()
events = events[events[:, 2] == 1][:2]
raw.load_data().pick_channels([raw.ch_names[pick] for pick in picks[::30]])
raw.info.normalize_proj()
del picks
sfreq_new = raw.info['sfreq'] / decim
with raw.info._unlock():
    raw.info['lowpass'] = sfreq_new / 12.  # suppress aliasing warnings
pytest.raises(ValueError, epochs.decimate, -1)
pytest.raises(ValueError, epochs.decimate, 2, offset=-1)
pytest.raises(ValueError, epochs.decimate, 2, offset=2)
for this_offset in range(decim):
    epochs = Epochs(raw, events, event_id,
                    tmin=-this_offset / raw.info['sfreq'], tmax=tmax,
                    baseline=None)
    idx_offsets = np.arange(decim) + this_offset
    for offset, idx_offset in zip(np.arange(decim), idx_offsets):
        expected_times = epochs.times[idx_offset::decim]
        expected_data = epochs.get_data()[:, :, idx_offset::decim]
        must_have = offset / float(epochs.info['sfreq'])
        assert (np.isclose(must_have, expected_times).any())
        ep_decim = epochs.copy().decimate(decim, offset)
        assert (np.isclose(must_have, ep_decim.times).any())
        assert_allclose(ep_decim.times, expected_times)
        assert_allclose(ep_decim.get_data(), expected_data)
        assert_equal(ep_decim.info['sfreq'], sfreq_new)

# More complicated cases
epochs = Epochs(raw, events, event_id, tmin, tmax)
expected_data = epochs.get_data()[:, :, ::decim]
expected_times = epochs.times[::decim]
for preload in (True, False):
    # at init
    epochs = Epochs(raw, events, event_id, tmin, tmax, decim=decim,
                    preload=preload)
    assert_allclose(epochs.get_data(), expected_data)
    assert_allclose(epochs.get_data(), expected_data)
    assert_equal(epochs.info['sfreq'], sfreq_new)
    assert_array_equal(epochs.times, expected_times)

    # split between init and afterward
    epochs = Epochs(raw, events, event_id, tmin, tmax, decim=dec_1,
                    preload=preload).decimate(dec_2)
    assert_allclose(epochs.get_data(), expected_data)
    assert_allclose(epochs.get_data(), expected_data)
    assert_equal(epochs.info['sfreq'], sfreq_new)
    assert_array_equal(epochs.times, expected_times)
    epochs = Epochs(raw, events, event_id, tmin, tmax, decim=dec_2,
                    preload=preload).decimate(dec_1)
    assert_allclose(epochs.get_data(), expected_data)
    assert_allclose(epochs.get_data(), expected_data)
    assert_equal(epochs.info['sfreq'], sfreq_new)
    assert_array_equal(epochs.times, expected_times)

    # split between init and afterward, with preload in between
    epochs = Epochs(raw, events, event_id, tmin, tmax, decim=dec_1,
                    preload=preload)
    epochs.load_data()
    epochs = epochs.decimate(dec_2)
    assert_allclose(epochs.get_data(), expected_data)
    assert_allclose(epochs.get_data(), expected_data)
    assert_equal(epochs.info['sfreq'], sfreq_new)
    assert_array_equal(epochs.times, expected_times)
    epochs = Epochs(raw, events, event_id, tmin, tmax, decim=dec_2,
                    preload=preload)
    epochs.load_data()
    epochs = epochs.decimate(dec_1)
    assert_allclose(epochs.get_data(), expected_data)
    assert_allclose(epochs.get_data(), expected_data)
    assert_equal(epochs.info['sfreq'], sfreq_new)
    assert_array_equal(epochs.times, expected_times)

    # decimate afterward
    epochs = Epochs(raw, events, event_id, tmin, tmax,
                    preload=preload).decimate(decim)
    assert_allclose(epochs.get_data(), expected_data)
    assert_allclose(epochs.get_data(), expected_data)
    assert_equal(epochs.info['sfreq'], sfreq_new)
    assert_array_equal(epochs.times, expected_times)

    # decimate afterward, with preload in between
    epochs = Epochs(raw, events, event_id, tmin, tmax, preload=preload)
    epochs.load_data()
    epochs.decimate(decim)
    assert_allclose(epochs.get_data(), expected_data)
    assert_allclose(epochs.get_data(), expected_data)
    assert_equal(epochs.info['sfreq'], sfreq_new)
    assert_array_equal(epochs.times, expected_times)

    # test picks when getting data
    picks = [3, 4, 7]
    d1 = epochs.get_data(picks=picks)
    d2 = epochs.get_data()[:, picks]
    assert_array_equal(d1, d2)

I meant that I thought I had addressed it but you're free to unresolve, sorry if that offended you.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it's not plucked from nowhere. 4 == evoked.times.tolist().index(0) - decim and 5 is because the second one has offset=1

Copy link
Member

@drammock drammock Jul 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

    evoked = EvokedArray(data, info, tmin=-1)
    zero_idx = evoked.times.tolist().index(0)
    evoked_dec = evoked.copy().decimate(decim)
    evoked_dec_2 = evoked.copy().decimate(decim, offset=1)
    evoked_dec_3 = evoked.decimate(dec_1).decimate(dec_2)
    start_samp = zero_idx - decim
    assert_array_equal(evoked_dec.data, data[:, start_samp::decim])
    # this has +1 because offset=1 when decimating ↓↓↓↓↓↓↓↓↓↓↓↓↓↓
    assert_array_equal(evoked_dec_2.data, data[:, (start_samp + 1)::decim])

you can change n_times on line 75 to whatever you like and this will still pass

mne/time_frequency/tests/test_tfr.py Show resolved Hide resolved
mne/utils/mixin.py Outdated Show resolved Hide resolved
mne/utils/mixin.py Outdated Show resolved Hide resolved
mne/utils/mixin.py Show resolved Hide resolved
@alexrockhill
Copy link
Contributor Author

That looks pretty sensible to me for the multiple return type https://output.circle-artifacts.com/output/job/76871b03-098a-4353-9ac0-402f230ccd3b/artifacts/0/dev/generated/mne.Evoked.html?highlight=evoked#mne.Evoked.decimate, definitely an improvement over main.

I'm not sure about the glossary entry for "time-series object". Try git grep "inst :" and you'll see lots of examples of how we handle this situation elsewhere. We're not perfectly consistent in how we do it, but let's not make that worse by adding a new variant.

I looked but I think this is a new case for mixins where you don't have multiple input types accepted (because it's a method) and ideally you would just put the correct output type as the return but there's no way to know the class from the method docstring so I think this is the most elegant solution rather than having a huge list or for me even more confusing evoked methods that say they return epochs.

@alexrockhill
Copy link
Contributor Author

alexrockhill commented Jul 22, 2022

Ok, I touched the failing examples and everything is green, should be good to go @drammock!

Unless you don't like the glossary entry but I looked pretty hard to find out how to get the method from the first inherited parent class and I don't think you can in python3 and so this is the best thing I could come up with which I think is better than what's on main (evoked returns epochs). The other option is a long list of return types but I think that's misleading as well because it will only return one type; the object that's method it is. I do really wish there were a way to infer the class from the method but it just doesn't appear to be possible, and the docstring is made only once for the parent class and then copied to the inherited classes so that seems to be insurmountable.

@alexrockhill
Copy link
Contributor Author

alexrockhill commented Jul 25, 2022

Ok to merge @drammock?

Thanks for writing the more comprehensive tests!

doc/conf.py Outdated
@@ -237,6 +237,7 @@
'Transform': 'mne.transforms.Transform',
'Coregistration': 'mne.coreg.Coregistration',
'Figure3D': 'mne.viz.Figure3D',
'MNE-object': ':term:`MNE time-series object`',
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd rather not do this in this PR. In most places we just do instance of Raw, Epochs, or Evoked. It's not optimal and in some cases not even complete (e.g., can be AverageTFR, STC, etc. instead) but it's common and people are used to it. If we want to go this route it should be in a separate PR where we change a ton of these at once. But for now I'd just make this instance of Raw, Epochs, Evoked, AverageTFR, or SourceEstimate or so.

But even separate-PR aside, I'm not totally convinced we need to or want to create this new category just for the purposes of making our Returns section better. For example, to me the best solution would be to somehow magically infer the class name being documented but I don't think this is possible based on how autodoc works. To me the simplest solution is to do something like instance of original type or something more compact that has this same meaning to say "you get back the same instance mostly for chaining purposes".

We should discuss this separately in a new issue. @alexrockhill if you're motivated by this problem, could you try making a new issue with the points that people have made about this already, then we can continue the discussion there?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1 for instance of original type, it is a tricky corner case for autodoc. Sounds good, I'll remove.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

... but for now in this PR please just use a list and open an issue where we discuss separately what the correct verbiage/linking should be for this "original type" stuff would be

larsoner added 2 commits July 27, 2022 10:06
* upstream/main:
  FIX: ctf data can sometimes have no date time. (mne-tools#10957)
  replace check_version from _testing.py with the one in check.py (mne-tools#10958)
  ENH:  mne-tools#10542 Adds __eq__() to DigMontage class and updates tests. (mne-tools#10942)
  Use HiDPI icons with PyQt5/PySide2 (mne-tools#10956)
@larsoner larsoner enabled auto-merge (squash) July 27, 2022 14:14
@larsoner larsoner added this to the 1.1 milestone Jul 27, 2022
@larsoner larsoner disabled auto-merge July 27, 2022 15:09
@larsoner larsoner merged commit 23f37cd into mne-tools:main Jul 27, 2022
@larsoner
Copy link
Member

Thanks @alexrockhill !

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants