ENH: Speed up set_bipolar_reference

[marsipu]

Reference issue

Fixes #8718

What does this implement/fix?

This improves the performance of adding a larger number (>25) of bipolar channels as suggested by @jasmainak.
EDIT: A gist to measure the speed can be found here

Additional information

I encountered some difficulties with how to best organize the creation of the new channels with the new method and handling Info while staying conservative. I think a large part of the performance improvement comes from moving add_channels outside the loop as done with the last commit 29f609d2. Before, even with the new matrix-multiplication-method the performance was much smaller with n>25.

At the moment, the new bipolar-channels are always appended (like when drop_refs=False) and don't take the place of the anode as in the previous implementation. If that is important I could add a reorder_channels afterwards with bipolar channels taking place of the anodes again.

As the bipolar-channels are not taking the same place as the anodes, I had to change the test a bit for now to pass.
There was also an additional attribute of the UpdateChannelsMixin called _projector besides picks which seemed to be outdated too with changed channels in Info.

Was it on purpose, that the bipolar reference of epochs was created with the mean of the anode-epochs?

ref_data = data[..., ref_from, :].mean(-2, keepdims=True)

In the new implementation it is created with the anode-signal of each single-epoch but could also be changed to the epochs mean.

ref_inst._data = np.asarray([multiplier.dot(ep) for ep in inst._data])

I changed the location of the bipolar-channel from zero to the location of the cathode as suggested by @jasmainak . I wondered, why every other 'chs'-attribute should still be taken from anode. Couldn't it be the anode-location instead or all the other info-attributes taken from the cathode?

What do you think about these changes?

[jasmainak]

couldn't you simply do:

ref_inst = inst.copy().pick_channels(cathode)
ref_inst.rename_channels(...)

it will read cleaner and maybe even faster

[marsipu]

I did this at first, but this fails when one channel occurs multiple times in anode or cathode, as pick_channels converts to a set. It would be cleaner and probably faster, but at the cost, that you can only use each channel once on each orientation

[marsipu]

There is a test for this in test_reference which fails then

[jasmainak]

ughgh ... I see. I would try to avoid using add_channels because that will probably make unnecessary copies of the data (for each new channel?). Can you use RawArray or EvokedArray or some such thing to construct the final instance in one go? I suspect it will be faster

[larsoner]

I would try to avoid using add_channels because that will probably make unnecessary copies of the data (for each new channel?). Can you use RawArray or EvokedArray or some such thing to construct the final instance in one go? I suspect it will be faster

Indeed iterating will be slower than a vectorized solution. I would figure out what dot product needs to be done, then do a single raw.drop_channels call, then do a single raw.add_channels call with the new data. The overhead from this should be much smaller

[jasmainak]

can you share some performance benchmarks? Say if I have 180 channels (so 179 anode-cathode pairs), how long does it take to compute bipolar with old and new implementation for a) raw b) epochs c) evoked

[jasmainak]

I would leave this for the end but there might be some performance benefits of using einsum for the epochs. I always get confused when using it though

[larsoner]

@ has much saner broadcasting behavior than np.dot (@ operates along the last two dims of both inputs, and np.dot... doesn't), this is probably the same as multiplier @ inst._data

[marsipu]

Thank you, @ works great. It doesn't seem to greatly improve performance though.

It appears to me, that the rereferencing itself doesn't take as much of the computation time as the overhead from creating and checking info in add_channels. You can try it out, if you replace l. 509-511 with:

    if isinstance(ref_instances, list):
        ref_instances = add_inst
    else:
        ref_instances.add_channels([add_inst], force_update_info=True)
ref_inst = ref_instances

[marsipu]

can you share some performance benchmarks? Say if I have 180 channels (so 179 anode-cathode pairs), how long does it take to compute bipolar with old and new implementation for a) raw b) epochs c) evoked

I uploaded a test-file for speed.
There is an error for epochs with this test-file using the old implementation which is now fixed.
Raw and Evokeds alone should demonstrate the drastic improve in performance too:

Before fix(-epochs):
Creation of 26 bipolar channels for Raw&Evoked took: 38.687 s

After fix (- epochs, - @-operator):
Creation of 26 bipolar channels for Raw&Evoked took: 1.801 s

After fix (+ epochs, - @-operator):
Creation of 26 bipolar channels for Raw&Evoked took: 3.241 s

After fix (+ epochs, + @-operator):
Creation of 26 bipolar channels for Raw&Evoked took: 3.218 s

I set the number to 26 because (at least on my laptop) the computation time seems to increase exponentially above ~20 (EDIT: and won't finish). I don't really understand why, I assume it is because of some kind of overhead introduced by the add_channels-code when called multiple times in a row.
Can you reproduce this behaviour?

[jasmainak]

Excellent, thanks for the benchmarks! You can edit your original description and add it there so it is not lost amongst the comments. And we can update it when you make more optimizations. 180 was the number of channels in my dataset and how I discovered the problem. You can see the dramatic improvement even with 26 channels :)

I think at this point, you can also try to use snakeviz or line profiler to see which lines are the slowest. Like you, I suspect add_channels though profiling will make it concrete

I think we should aim for sub-second speed, referencing should be instantaneous!

ps: it might be easier if you share a gist that works with sample data (than a custom file) so we can copy-paste instantly and try things if you need help

[marsipu]

A quick question:
I don't really understand, why electrode-information is taken from the anode, but location should be taken from the cathode.
Is this a convention?

[agramfort]

it should not be. the chs['loc'] should ideally contain the 2 locations with the 2nd location being the reference used.

…

[marsipu]

it should not be. the chs['loc'] should ideally contain the 2 locations with the 2nd location being the reference used.
…

Ok thank you, now all channel-info including location is copied from anode

[marsipu]

In this gist I uploaded the protocols of profiling before and after the most recent improvements. The major part of the time was indeed taken by the multiple call of add_channels and the creation of an instance for each new bipolar-channel.
Now, as suggested, only one instance for the new rereferenced channels is created which is then added once to the original instance. I rather kept the add_channels once, because I didn't want to miss anything related to merging data and info, even though the two instances should be very similar. Or do you think this should be improved too with a custom function for adding the referenced data?

Performance before recent improvements:

Creation of 26 bipolar channels took:
Raw: 1.642 s
Epochs: 2.073 s
Evoked: 0.803 s

Creation of 180 bipolar channels took:
Raw: 11.286 s
Epochs: 14.424 s
Evoked: 5.513 s

Performance after recent improvements:

Creation of 26 bipolar channels took:
Raw: 0.165 s
Epochs: 0.275 s
Evoked: 0.054 s

Creation of 180 bipolar channels took:
Raw: 0.358 s
Epochs: 0.633 s
Evoked: 0.055 s

[drammock]

@marsipu what is the reason for [ci skip] here? The code changes are substantial, they should be tested. Can you git commit --amend to remove the [ci skip] from your last commit, and then git push --force? Or alternatively, git commit --allow-empty -m "trigger CIs" and then git push.

[marsipu]

@marsipu what is the reason for [ci skip] here? The code changes are substantial, they should be tested. Can you git commit --amend to remove the [ci skip] from your last commit, and then git push --force? Or alternatively, git commit --allow-empty -m "trigger CIs" and then git push.

@drammock Oh I am sorry, I thought it was good practice to skip ci for minor changes to save resources for other PRs. But I see now, I probably should have had at least documentation ci again for cb2e3d0. There was ci for all major changes before. When is it allowed/advisable to skip any ci?
I will trigger ci again as you instructed.

[drammock]

When is it allowed/advisable to skip any ci?

It is generally speaking not avisable to skip CIs. The one exception that I'm willing to put in writing is this: if all you changed is documentation (tutorial, example, and/or docstring), then [skip azp][skip github] is acceptable.

The reason for this is that the CIs can catch problems unrelated to your changes, that cropped up due to new versions of our dependencies being released. The sooner we catch those, the sooner we can fix them before real users encounter them. Two such examples came up just yesterday: #9321 and pydata/pydata-sphinx-theme#395.

The other case where you might do [ci skip] is if your PR is a draft PR and you know that the commits you're pushing aren't yet enough to fix the bug / implement the new feature, and you don't want to waste energy by running the CIs on something you already know will fail.

[marsipu]

It is generally speaking not avisable to skip CIs. The one exception that I'm willing to put in writing is this: if all you changed is documentation (tutorial, example, and/or docstring), then [skip azp][skip github] is acceptable.

Thank you for clarification, I will respect that from now on.

[larsoner]

LGTM +1 for merge with or without my test suggestion. Thanks @marsipu !

[jasmainak]

Fantastic PR. Thanks @marsipu !

[marsipu]

There seems to be an error for test_report.py::test_scraper, which doesn't seem to be related to this PR, right?
I could maybe fix this with exist_ok=True for makedirs in the line throwing the error, but I am not sure if this defies the purpose of the test.

[drammock]

I've opened #9332 to deal with the failing scraper tests.

[larsoner]

Close-reopen cycle to restart all CIs

[larsoner]

Thanks @marsipu !