[MRG] ENH: Add weighted symbolic mutual information (wSMI) connectivity measure #307
Conversation
GioMarraffini
force-pushed
the
wsmi-feature
branch
from
d98af9b to
755fe3e
Compare
GioMarraffini
force-pushed
the
wsmi-feature
branch
2 times, most recently
from
907406d to
e9cdf80
Compare
GioMarraffini
changed the title
GioMarraffini
force-pushed
the
wsmi-feature
branch
from
e9cdf80 to
efc4c45
Compare
tsbinns
left a comment
tsbinns
left a comment
There was a problem hiding this comment.
Thanks for opening the PR! Below are a couple suggestions/comments.
Typically, new features like this would also include an example. This can be just a very short demonstration of using the function (e.g., https://mne.tools/mne-connectivity/dev/auto_examples/mne_inverse_envelope_correlation_volume.html) or it could be a more detailed tutorial with an explanation of the method (e.g., https://mne.tools/mne-connectivity/dev/auto_examples/dpli_wpli_pli.html).
I haven't checked the unit tests extensively yet, will revisit based on the discussions. From a first glance they look extensive and the regression tests are great to have.
|
Just added the changes to a new commit with the example in mne-connectivity/examples/wsmi_example.py. I dont know if this is correct please let me know |
GioMarraffini
force-pushed
the
wsmi-feature
branch
from
194f51a to
c15e2d2
Compare
|
@GioMarraffini @Laouen Thanks both for responding to my initial suggestion so quickly. The changes with the anti-aliasing parameter look very reasonable. |
@tsbinns Hi just wanted to ask if you had time to review it and if there were any news. We are happy to make any necessary changes. |
tsbinns
left a comment
tsbinns
left a comment
There was a problem hiding this comment.
Sorry for the delay. Altogether I think the changes were very good and I trust your judgement for adapting the anti-aliasing parameter. Here are a few more comments from me.
|
Rerunning the failing test; should not be a problem from your end. |
|
@tsbinns Thanks for the revisions and additional comments! I've just pushed the new updates. let me know if you spot anything else. |
tsbinns
left a comment
tsbinns
left a comment
There was a problem hiding this comment.
Thanks for the revisions! Some very minor points here. The main outstanding things are:
- Supporting arrays as data inputs from before.
- Triaging numba to be an optional dependency.
- Switching to catching memory errors rather than a user-specified memory limit (pending input from others).
It's really great work so far! Let me know if you need a hand bringing any of this over the line.
|
Something I also didn't consider before is the size of the testing data. The existing test data is 0.9 MB, but the wSMI data is 5.3 MB. Without the wSMI test data the source directory is only ~3 MB. |
|
The problem with generating test data online is that we have no ground thruth for that data. Currently, the ground truth is the value obtained by the original function. Maybe we can do some tests with estimated result patterns and if we see those partterns then assume is correct even do we don't compare with exact results. |
So the way we have handled this in the past when porting methods from other toolboxes is to have a regression test involving some arbitrary connectivity data of only a few channels/epochs and comparing the results from the original and our implementations. This way the data we need to save is very small.
This is what we have done for evaluating whether ported methods are working as intended. We simulate test data in the scripts with various connectivity profiles and validate the results based on whether the connectivity values are in a reasonable range. We don't have to check that they match the original implementation exactly since the regression test is handling that side of things.
|
|
Hi @GioMarraffini, do you think those latest commits in a8d8263 are ready for another round of reviews? |
Hi @tsbinns. Yes, Thanks! If you see something to change please let me know |
There was a problem hiding this comment.
Hi @GioMarraffini, overall the changes were great and adressed the previous comments. There are just a few things regarding the array inputs.
Also, we don't get notifications for commits, so just add a comment or ping us once you think the code is ready for review.
|
Hi @tsbinns , I've addressed all the feedback from your last review. The code is ready for another check. |
tsbinns
left a comment
tsbinns
left a comment
There was a problem hiding this comment.
@GioMarraffini Just some very minor things again.
|
Hi @tsbinns , I've addressed your latest feedback. The implementation now follows the pattern you suggested. Ready for review! |
Looks good! |
drammock
left a comment
drammock
left a comment
There was a problem hiding this comment.
will resume review next week.
mne_connectivity/tests/test_wsmi.py
Outdated
| ) | ||
|
|
||
| # Results should be identical | ||
| assert_allclose(conn1.get_data(), conn2.get_data()) |
There was a problem hiding this comment.
if they should be truly identical, use assert_array_equal
mne_connectivity/tests/test_wsmi.py
Outdated
| assert conn.method == "wSMI" | ||
|
|
||
|
|
||
| def test_wsmi_all_channels_as_bad(): |
There was a problem hiding this comment.
you're already testing what happens if 1 channel is bad, and since the behavior is apparently to ignore bad channels, we shouldn't bother also testing when 2, 3, 4, or all channels are bad.
mne_connectivity/tests/test_wsmi.py
Outdated
| assert conn.method == "wSMI" | ||
| assert conn.n_nodes == n_channels | ||
| assert np.all(np.isfinite(conn.get_data())) |
There was a problem hiding this comment.
A test that parametrizes over kernel and tau values should test that the different kernel or tau values have different effects. There are at least 2 ways to do that:
- use non-random channel data, compute from first principles what the conn values should be for different kernel/tau values, pass them in as part of the parametrization, and assert that you get what you expect
- don't parametrize, run the wsmi a few times (with different kernel/tau values) within the test function, and assert that the resulting
conn.get_data()arrays are actually different.
mne_connectivity/tests/test_wsmi.py
Outdated
| # Values should be different (wSMI uses weights, SMI doesn't) | ||
| # They shouldn't be exactly equal due to the weighting | ||
| # The methods differ in their weighting, but with random data | ||
| # they might produce similar results (this is expected behavior) | ||
| # So we just check that both produce finite values | ||
| pass |
There was a problem hiding this comment.
I think the right lesson to draw from this is that we shouldn't be using random data in tests like this one.
mne_connectivity/tests/test_wsmi.py
Outdated
| # ============================================================================= | ||
| # Ground Truth Validation Tests (New Test Data System) | ||
| # ============================================================================= |
There was a problem hiding this comment.
ah, I wish I'd known this was coming! I'm going to stop reviewing for now because I'm out of time today. I haven't yet reviewed the tests that follow, so take this with a grain of salt, but most of the preceding tests can probably be removed or consolidated. E.g, there could maybe be 1 test using fake/random data (could even be np.zeros) that checks all the expected errors/warnings given different degenerate input values; maybe 1 test to verify a bunch of basic shape/object properties, and then everything else should probably be these "ground truth" type of tests (if done correctly, they ought to hit most/all of the variations you hit above, e.g., different channel types, weighting, antialias, averaging, time-windowing, or indexing).
There was a problem hiding this comment.
Hi @drammock ! thanks for all the feedback! just one question should I wait until these following reviews on this file or should I start deleting previous tests? I already made the necesary adjustment to the example file. Please let me know if it needs further modifications. Thanks!
There was a problem hiding this comment.
Yes feel free to start simplifying the tests now. I'll try to take another look on Friday
GioMarraffini
force-pushed
the
wsmi-feature
branch
from
e962b81 to
6936b89
Compare
GioMarraffini
force-pushed
the
wsmi-feature
branch
from
cd77750 to
01cbd63
Compare
|
Hi @GioMarraffini, I am in the middle of double checking the bad channel handling behaviour, which I will (really) try and finish this week, but I would push any needed changes so that would not require more time from you. |
Thanks @tsbinns . Sorry for the rush, but since @GioMarraffini internship at the lab will finish soon, it would be a shame not to be able to merge this PR. Thanks again for all your time invested on this PR! |
drammock
left a comment
drammock
left a comment
There was a problem hiding this comment.
Hi @GioMarraffini and @Laouen,
Apologies for the very long delay in reviewing. @tsbinns and I just synced up about this PR, and have concluded that there's a fair bit of work left to get this PR into a state where we're comfortable merging it. As far as we can tell the implementation is fine, and you've done a thorough job generating the example and the tests. However, we feel the tests are too long/numerous and also not stringent enough. There are also some testing "tricks" we'd like to use (parametrizing test inputs, possibly fixtures for the simulated data) that would make this more maintainable in the long-term.
I understand that your internship ends soon and you would like to have this all-the-way done before then, but we simply cannot guarantee that timeline; reviewing large PRs takes a long time, and we have to be judicious about the added maintenance burden we're taking on here. If you want to keep working on it (now or after your internship) you can of course do so, but FYI our plan is for Thomas and I to push some commits in the coming week or so, to get things to a mergable state, so you can rest assured that the work will get merged in eventually (hopefully fairly soon!) even if you want/need to stop working on it now.
(NB: inline comments below are somewhat terse, as they're meant more as notes-to-self than as feedback directed at the authors)
examples/wsmi.py
Outdated
|
|
||
| for epoch in range(n_epochs): | ||
| # Set random seed for reproducibility within epoch variation | ||
| rng = np.random.default_rng(42 + epoch) |
| # 4. Since channel 1 is identical to channel 0, its coupling to channel 2 | ||
| # should be the same as channel 0's coupling to channel 2 | ||
| assert abs(coupled_wsmi - identical_coupled_wsmi) < 1e-10, ( | ||
| "wSMI should be identical for identical source channels" | ||
| ) |
There was a problem hiding this comment.
shouldn't it be possible to assert true equality here? the signals in ch1 and ch2 are identical... is WSMI somehow non-deterministic?
| # 6. Test reproducibility across epochs | ||
| assert np.all(np.std(conn_data, axis=0) < 0.5), ( | ||
| "wSMI should be stable across epochs" | ||
| ) |
There was a problem hiding this comment.
TODO check stringency is appropriate
mne_connectivity/tests/test_wsmi.py
Outdated
| # Create test data with clear nonlinear coupling | ||
| data = np.zeros((n_epochs, 3, n_times)) | ||
|
|
||
| for epoch in range(n_epochs): |
There was a problem hiding this comment.
should be able to simulate without a for-loop?
mne_connectivity/tests/test_wsmi.py
Outdated
| # Test tau=1 vs tau=2 | ||
| conn_tau1 = wsmi(epochs, kernel=3, tau=1) | ||
| conn_tau2 = wsmi(epochs, kernel=3, tau=2) | ||
|
|
||
| data_tau1 = np.mean(conn_tau1.get_data(), axis=0) | ||
| data_tau2 = np.mean(conn_tau2.get_data(), axis=0) | ||
|
|
||
| # Connection indices: (0,1)=0, (0,2)=1, (1,2)=2 | ||
| nonlinear_tau1 = data_tau1[0] # base-nonlinear with tau=1 | ||
| independent_tau1 = data_tau1[1] # base-independent with tau=1 | ||
| nonlinear_tau2 = data_tau2[0] # base-nonlinear with tau=2 | ||
| independent_tau2 = data_tau2[1] # base-independent with tau=2 | ||
|
|
||
| # Test that tau=2 shows better discrimination than tau=1 | ||
| if independent_tau1 > 0 and independent_tau2 > 0: | ||
| ratio_tau1 = nonlinear_tau1 / independent_tau1 | ||
| ratio_tau2 = nonlinear_tau2 / independent_tau2 | ||
|
|
||
| # tau=2 should show better discrimination (higher ratio) | ||
| assert ratio_tau2 > ratio_tau1, ( | ||
| f"tau=2 should show better nonlinear discrimination: " | ||
| f"tau=1 ratio={ratio_tau1:.2f}, tau=2 ratio={ratio_tau2:.2f}" | ||
| ) | ||
|
|
||
| # tau=2 should show substantial improvement (at least 2x better) | ||
| assert ratio_tau2 > 2 * ratio_tau1, ( | ||
| f"tau=2 should show substantial improvement over tau=1: " | ||
| f"tau=2 ratio should be > 2x tau=1 ratio" | ||
| f"tau=2 ratio / tau=1 ratio: {ratio_tau2 / ratio_tau1:.2f}" | ||
| ) |
mne_connectivity/tests/test_wsmi.py
Outdated
| assert np.all(np.isfinite(data_matrix)) | ||
|
|
||
| # Test with different channel types (suppress unit change warning) | ||
| import warnings |
There was a problem hiding this comment.
import should be at top of file... but don't do it that way either. Use with pytest.warns(...)
mne_connectivity/tests/test_wsmi.py
Outdated
| with warnings.catch_warnings(): | ||
| warnings.simplefilter("ignore", RuntimeWarning) | ||
| epochs.set_channel_types({"Fz": "eeg", "Cz": "eeg", "Pz": "mag", "Oz": "mag"}) | ||
| conn_mixed = wsmi(epochs, kernel=3, tau=1) | ||
| assert conn_mixed.n_nodes == 4 | ||
| assert np.all(np.isfinite(conn_mixed.get_data())) |
There was a problem hiding this comment.
move this check into test_wsmi_input_validation_and_errors
mne_connectivity/tests/test_wsmi.py
Outdated
| # Check basic properties | ||
| assert conn.method == "wSMI" | ||
| assert conn.n_nodes == 4 | ||
| assert conn.n_epochs_used == 2 | ||
|
|
||
| # Check data shape and validity | ||
| data_matrix = conn.get_data() | ||
| expected_connections = 4 * 3 // 2 # 4 choose 2 = 6 connections | ||
| assert data_matrix.shape == (2, expected_connections) | ||
| assert np.all(np.isfinite(data_matrix)) |
There was a problem hiding this comment.
I think all this is basic enough that we can skip checking it, or move it into the input_validation_and_errors test (possibly renamed to input_and_output_validation or so)
mne_connectivity/tests/test_wsmi.py
Outdated
| assert np.all(np.isfinite(conn_mixed.get_data())) | ||
|
|
||
|
|
||
| def test_wsmi_averaging_and_indices(): |
There was a problem hiding this comment.
I think this whole test can go away, except for the "invalid indices" error checks.
|
@GioMarraffini One comment looking at the example, what would the interpretation be of negative values? There are some occurring in the simulated data. 
|
|
To summarise those changes pushed in e93edca:
|
Individual terms can be negative: For a specific symbol pair (i,j), when the observed joint probability P(X=i, Y=j) is less than what you'd expect under independence P(X=i) * P(Y=j), that term contributes negatively to MI: Standard MI is always non-negative because when you sum over ALL symbol pairs, positive and negative contributions balance out (by Jensen's inequality) to give a non-negative result. However, wSMI breaks this balance: The weighting matrix (from _get_weights_matrix) selectively zeros out certain symbol pairs: What Negative wSMI Values Mean for Connectivity TL;DR: In the original sense, wSMI isn’t a signed “positive vs. negative connectivity” measure. Standard mutual information is ≥ 0, but wSMI applies a weighting that drops “trivial” symbol pairings (identical or opposite-sign) to reduce volume-conduction effects. That weighting breaks the usual KL-divergence guarantee, so small negative values can appear—they don’t mean “anticoupling”; they usually mean “no information sharing above chance (and/or estimator noise). |
|
@GioMarraffini Thanks for the detailed answer! That makes sense. I think it would be good to include a short note about that in the example, perhaps after the results are shown in the Visualizing wSMI Connectivity Results section. Could be something like:
What do you think? |
Great, I think that would help |
|
Hello @tsbinns here I answer regarding your tests comment:
I think we are aligned and your description of regresion tests is what we have currently impemented, with the unique distinction that we've used the whole original data instead of a minimum required tests set as you have in the data you point out.
perfect, we can do this kind of tests as well. We will implement a tests in this direction and push a new commit. |
Thanks a lot for your detailed feedback and for taking such careful attention with this PR — it’s great to see how much care goes into maintaining the library’s standards. We completely agree that ensuring quality is the top priority, even if it takes more time. That said, since @GioMarraffini’s internship is ending soon, we were wondering if it might be possible to organize a short “PR party” or Zoom call in the coming days. This could help us go through the code together and speed up the review process a bit — without, of course, compromising the quality and rigor of the review. |
|
@Laouen @GioMarraffini We can arrange a call to go through the PR together. The best time would be during MNE's Office Hours held every other Friday. @drammock would be there, and I could also join. The next session is 7th Nov. |
@tsbinns @drammock This friday 7th at 17h00 is good for us as we see that is your Office Hours. we'are looking forward to see you this friday! Best |
|
Hey @Laouen @GioMarraffini, we are around at the moment in the office hours. Feel free to drop by anytime until the end of the hour. |
Hi @tsbinns sorry for not being able to connect in your last office hours. Sadly @GioMarraffini is currently on a conference so we can't yet have a time. Can you give us the calendar with the next office hours so we organize ourself to make a slot and finally do that call ? |
| # Plot a sample of the data | ||
| fig = epochs.plot( | ||
| n_epochs=1, | ||
| scalings="auto", | ||
| show_scrollbars=False, | ||
| title="Synthetic EEG Data with Different Connectivity Patterns", | ||
| ) |
There was a problem hiding this comment.
Plot not rendering in docs, will need to resolve before merge.
Hi @Laouen, the next office hours isn't until the 19th (calendar link here: https://mne.discourse.group/t/mne-office-hours-on-discord/7439). |
@tsbinns Perfect, thanks for the calendar link. @GioMarraffini and I are going to connect in the next office hours (12/19) to the discor so we can have the meeting with you. Looking forward to discuss this PR with you! |
5 participants
Thanks for contributing. If this is your first time,
make sure to read contributing.md
PR Description
Add wsmi() function to compute weighted symbolic mutual information, a connectivity measure for studying markers of consciousness and conscious perception based on symbolic dynamics.
The implementation includes: - Symbolic transformation using ordinal patterns with optimizations - Mutual information computation between symbolic sequences using Numba JIT - Weighting based on pattern distance for enhanced sensitivity - Memory checks for large kernel sizes to prevent excessive memory usage - Comprehensive test suite with ground truth validation against reference data - Integration with MNE-Connectivity's EpochTemporalConnectivity class - Support for Current Source Density (CSD) preprocessing for EEG data - Configurable low-pass filtering with automatic frequency selection
The measure quantifies non-linear statistical dependencies between time series and is particularly useful for studying consciousness-related brain dynamics.
This PR closes <#301> (Issue number)
Merge checklist
Maintainer, please confirm the following before merging: