-
Notifications
You must be signed in to change notification settings - Fork 1.4k
[MRG] PCA flip for volumetric source estimates #13092
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
base: main
Are you sure you want to change the base?
Conversation
Extension of PCA flip in volumetric setting. Underlying logic is pretty much identical to cortical case with respect to flip vector creation and PCA itself Added a check in PCA flip in case the flip is None or number of vertices is below 2, in which case PCA will return a trivial estimate of the signal Complete with its own unit test case
I belive this PR is ready for review :) (although I do not know why pip pre is failing on Windows) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Found some cruft, can you take a look at the diff
on GitHub and make sure I found all of it?
Also with enhancements it's nice to see the effect in some example. Maybe you could modify this:
https://mne.tools/1.8/auto_examples/inverse/compute_mne_inverse_epochs_in_label.html
which comes from examples/inverse/compute_mne_inverse_epochs_in_label.py
could be modified to add a volume example to show its flip modes as well? Or maybe some other example like tutorials/inverse/60_visualize_stc.py
or examples/inverse/label_source_activations.py
?
mne/source_estimate.py
Outdated
# logger.info("Selected mode: " + mode) | ||
# print("Entering _prepare_label_extraction") | ||
# print("Selected mode: " + mode) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Cruft, should be removed. But really to help next time:
# logger.info("Selected mode: " + mode) | |
# print("Entering _prepare_label_extraction") | |
# print("Selected mode: " + mode) | |
logger.debug(f"Selected mode: {mode}") |
then if you run with verbose='debug'
you will see the statement
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed with logger.debug, thanks!
mne/source_estimate.py
Outdated
@@ -3445,6 +3456,7 @@ def _prepare_label_extraction(stc, labels, src, mode, allow_empty, use_sparse): | |||
|
|||
bad_labels = list() | |||
for li, label in enumerate(labels): | |||
# print("Mode: " + mode + " li: " + str(li) + " label: " + str(label)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
# print("Mode: " + mode + " li: " + str(li) + " label: " + str(label)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed
mne/source_estimate.py
Outdated
# print(f"src: {src[:2]}") | ||
# print(f"len(src): {len(src[:2])}") | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
# print(f"src: {src[:2]}") | |
# print(f"len(src): {len(src[:2])}") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fixed
@@ -1460,22 +1460,47 @@ def label_sign_flip(label, src): | |||
flip : array | |||
Sign flip vector (contains 1 or -1). | |||
""" | |||
if len(src) != 2: | |||
raise ValueError("Only source spaces with 2 hemisphers are accepted") | |||
if len(src) > 2 or len(src) == 0: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A better / more modern check would be something like:
_validate_type(src, SourceSpaces, "src")
_check_option("source space kind", src.kind, ("volume", "surface"))
if src.kind == "volume" and len(src) != 1:
raise ValueError("Only single-segment volumes, are supported, got labelized volume source space")
And incidentally I think eventually we could add support for segmented volume source spaces, as well as mixed source spaces (once surface + volume are fully supported, mixed isn't too bad after that). But probably not as part of this PR!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Probably! To support this, the (clarified) code in label_sign_flip with the hemis dictionary might be recycled:
hemis = {}
# Build hemisphere info dictionary
if label.hemi == "both":
hemis["lh"] = {"id": 0, "vertno": src[0]["vertno"]}
hemis["rh"] = {"id": 1, "vertno": src[1]["vertno"]}
elif label.hemi in ("lh", "rh"):
hemis[label.hemi] = {"id": 0, "vertno": src[0]["vertno"]}
else:
raise Exception(f'Unknown hemisphere type "{label.hemi}"')
mne/label.py
Outdated
if isbi_hemi: | ||
lh_id = 0 | ||
rh_id = 1 | ||
lh_vertno = src[0]["vertno"] | ||
rh_vertno = src[1]["vertno"] | ||
elif label.hemi == "lh": | ||
lh_vertno = src[0]["vertno"] | ||
elif label.hemi == "rh": | ||
rh_id = 0 | ||
rh_vertno = src[0]["vertno"] | ||
else: | ||
raise Exception(f'Unknown hemisphere type "{label.hemi}"') |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't totally follow this logic. It's a surface source space if and only if it's bihemi, and it's volume source space if and only if it's a single-element list. I think it would be cleaner to triage based on source space type probably...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've clarified it to be clearer:
hemis = {}
# Build hemisphere info dictionary
if label.hemi == "both":
hemis["lh"] = {"id": 0, "vertno": src[0]["vertno"]}
hemis["rh"] = {"id": 1, "vertno": src[1]["vertno"]}
elif label.hemi in ("lh", "rh"):
hemis[label.hemi] = {"id": 0, "vertno": src[0]["vertno"]}
else:
raise Exception(f'Unknown hemisphere type "{label.hemi}"')
# get source orientations
ori = list()
for hemi, hemi_infos in hemis.items():
# When the label is lh or rh, get vertices directly
if label.hemi == hemi:
vertices = label.vertices
# In the case where label is "both", get label.hemi.vertices
# (so either label.lh.vertices or label.rh.vertices)
else:
vertices = getattr(label, hemi).vertices
vertno_sel = np.intersect1d(hemi_infos["vertno"], vertices)
ori.append(src[hemi_infos["id"]["nn"][vertno_sel]])
if len(ori) == 0:
raise Exception(f'Unknown hemisphere type "{label.hemi}"')
I believe the PR is ready for a second round :) |
some test failures look related: could you run these tests and debug locally? test_label_sign_flip, test_extract_label_time_course_volume_pca_flip, test_extract_label_time_course, test_extract_label_time_course_equiv, test_label_extraction_subject |
All offending tests should now correctly be passing (locally no more issue on these tests) |
], | ||
) | ||
def test_extract_label_time_course_volume_pca_flip( | ||
src_volume_labels, label_type, mri_res, test_label, cf, call |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We always want tests to be fairly complete, but competing against that we also want them to be as quick as possible. Our test suite takes over an hour to run already! 😢
Do we need to iterate over cf
here for example? If we already test that functionality elsewhere and the pca_flip
functionality here is independent of that, let's not. Same thing with mri_resolution
, this seems like something that is unrelated and tested elsewhere most likely.
Maybe the pca_flip
option could just be added somewhere else (maybe as one parameter in an already parametrized function) to really trim it down?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
On this one, I'm on the fence - integrating pca_flip as an option in another test is certainly an option, much like "mean" for example, but at the same time we lose clarity in the testing if we over parameterize no? In particular, it becomes tedious to know what is actually being tested and which cases might or might not be covered...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In this case I think you're over-testing quite a bit. Your code change is only relevant to some minor parts of the code, not coordinate frame, mri resolution, etc. So I think it makes more sense to trim it down (a lot!) or integrate elsewhere
What does this implement/fix?
This PR implements PCA flip in volumetric source estimates.
To do so, it relies on using the already existing _pca_flip method. _pca_flip has been slighly modified to check if at least two vertices are available within a label. If that is not the case, it returns trivially the original signal, as PCA is ill-defined in this case.
Additional information
Feel free to have a look and - of course - double-check for correctness!