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

ENH: Support TD data #11064

Draft
wants to merge 40 commits into
base: main
Choose a base branch
from
Draft

Conversation

larsoner
Copy link
Member

@larsoner larsoner commented Aug 21, 2022

Closes #9661

@rob-luke how about this plan:

For now I'm inclined to make it so if you get a SNIRF_PROCESSED data type that we don't know how to handle, we should maybe raise an error in this case...? For now I've removed it as a coil_type because this _PROCESSED really seems to mean something else specific, and we should perhaps add those one by one as we find use cases for them.

I have some collaborators who have seen Kernel help files telling them to use this very out of date branch, so it would be nice to get this in!

@larsoner larsoner added the fNIRS label Aug 23, 2022
@rob-luke
Copy link
Member

Thanks for jumping in to finish this off @larsoner

For now I'm inclined to make it so if you get a SNIRF_PROCESSED data type that we don't know how to handle, we should maybe raise an error in this case...?

Agreed

Ping me when you want a review, it doesn't need to be a final review, I am happy to take a look early if that helps.

@larsoner
Copy link
Member Author

Feel free to look now @rob-luke to see if you're happy with the overall design, the rest of the work is "just" some remaining details (adding unit tests, playing with the data a little bit to set scale factors, etc.) that I think we can rely on CIs to check, or that we can iterate on in follow-ups (e.g., tweaking plot scale factors)

Copy link
Member

@rob-luke rob-luke left a comment

Choose a reason for hiding this comment

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

Looking good, it will be great to have support for TD in MNE! I will review again once there are tests.

mne/defaults.py Outdated Show resolved Hide resolved
larsoner added 4 commits March 1, 2023 12:46
* upstream/main: (264 commits)
  BUG: Fix deprecated API usage in example (mne-tools#11512)
  Deprecate 'kind' and 'path' in favor of 'fname' in the layout reader (mne-tools#11500)
  EGI/MFF events outside EEG recording should not break the code (mne-tools#11505)
  fixed annotations error on export (mne-tools#11435)
  DOC: Update installer links [skip azp] [skip actions] [skip cirrus] (mne-tools#11506)
  BUG: updates for MPL 3.7 compatibility (mne-tools#11409)
  Fix docstrings by replacing str with path-like and fix double backticks for formatting (mne-tools#11499)
  Use pathlib.Path instead of os.path to handle files and folders [circle deploy] (mne-tools#11473)
  MAINT: Fix Circle [circle deploy] (mne-tools#11497)
  MAINT: Use mamba in CIs (mne-tools#11471)
  Updating documentation to clarify full vs half-bandwidth and defaults in time_frequency.multitaper.py (mne-tools#11479)
  Fix typo in tutorial (mne-tools#11498)
  Typo fix and added colons before code (mne-tools#11496)
  [MRG] ENH/DOC: demo custom spectrum creation (mne-tools#11493)
  Accept only left-clicks for adding annotations (mne-tools#11491)
  [BUG, MRG] Fix pial surface loading, logging in locate_ieeg (mne-tools#11489)
  [ENH] Added unit_role to add non-breaking space between magnitude  and units (mne-tools#11469)
  MAINT: Fix CircleCI build (mne-tools#11488)
  [DOC] Updated decoding.SSD documentation and internal variable naming (mne-tools#11475)
  Typo fix (mne-tools#11485)
  ...
larsoner and others added 3 commits October 28, 2024 13:25
Co-authored-by: Robert Luke <[email protected]>
Co-authored-by: John Griffiths <[email protected]>
Co-authored-by: Julien Dubois <[email protected]>
@@ -137,6 +137,7 @@
.. _Joan Massich: https://github.com/massich
.. _Johann Benerradi: https://github.com/HanBnrd
.. _Johannes Niediek: https://github.com/jniediek
.. _John Griffiths: https://www.grifflab.com
Copy link
Member Author

Choose a reason for hiding this comment

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

@JohnGriffiths @julien-dubois-k @Zahra-M-Aghajan let me know if you want a different URL to link from your name in the changelog

Choose a reason for hiding this comment

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

grifflab.com is perfect, thanks

Copy link
Member Author

@larsoner larsoner left a comment

Choose a reason for hiding this comment

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

Okay code has been updated. I suspect the scaling could be wrong as I mention in a comment above ☝️ below 👇 but in any case people can test it with some data!

It's possible/likely that in addition to the possible scaling bug, we'll want to adjust the default plot scaling limits.

It would also be great to have an example added (maybe using mne-misc-data) showing that this data can at least be read!

Comment on lines +446 to +447
fnirs_time_delays[bin_idx - 1]
* fnirs_time_delay_widths[bin_idx - 1]
Copy link
Member Author

Choose a reason for hiding this comment

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

I modified this to index into fnirs_time_delay_widths when it didn't before. It didn't error in either case, and it seems like one of them should have! So we still need to add a suitable SNIRF test file to mne-testing-data @JohnGriffiths (or others) if you have a tiny one you can share

Copy link
Member Author

Choose a reason for hiding this comment

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

@julien-dubois-k this one is also for you ☝️

Choose a reason for hiding this comment

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

test file shared here, hopefully it helps with this case

data = raw.get_data("fnirs_td_moments_amplitude")
assert data.shape == raw._data.shape
norm = np.nanmedian(np.linalg.norm(data, axis=-1))
assert 1e5 < norm < 1e6 # TODO: 429256, is this reasonable Molars!??
Copy link
Member Author

@larsoner larsoner Oct 28, 2024

Choose a reason for hiding this comment

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

This should be in units of M (mol). Looking at the values in the test file:

>>> import h5py
>>> import numpy as np
>>> x = h5py.File("td_moments.snirf", "r")
>>> list(x["nirs/data1/measurementList1"])  # no dataUnit, so code assumes M not uM!
['dataType', 'dataTypeIndex', 'dataTypeLabel', 'detectorIndex', 'sourceIndex', 'wavelengthIndex']
>>> np.array(x["nirs/data1/measurementList1/dataTypeLabel"]).item()
b'Time Domain - Moments - Amplitude'
>> np.nanmedian(np.linalg.norm(np.array(x["nirs/data1/dataTimeSeries"]), axis=0))
np.float64(429256.98709975014)

should the code assume that if no dataUnit is present, the units are actually uM?

Copy link
Member Author

Choose a reason for hiding this comment

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

@julien-dubois-k this one is for you ☝️

Copy link

@julien-dubois-k julien-dubois-k Nov 7, 2024

Choose a reason for hiding this comment

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

b'Time Domain - Moments - Amplitude' so this is moments data.
The units will be:

  • unitless for the "sum" moment (it's a photon count)
  • some multiple of s for the "mean" moment (it's a time of flight)
  • some multiple of s**2 for the "var" moment (it's a variance of the time of flight)

The SNIRF spec doesn't require dataUnit to be populated (here). In cedalion it looks like they just treat unspecified dataUnit as unitless or a.u. (here). I think vendors who generate SNIRF files should populate the dataUnit field always if there is a unit. Kernel didn't originally (hence the missing dataUnit in the current test files), but does now. I'll provide new test files shortly.

Choose a reason for hiding this comment

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

test files shared here

* upstream/main:
  eegbci api: allow downloading multiple subjects (mne-tools#12918)
  DOC: Document Linux desktop workaround (mne-tools#12900)
  Allow exporting edf where a channel contains only constant values (mne-tools#12911)
  Autogenerate environment.yml file from pyproject.toml (mne-tools#12914)
@larsoner larsoner marked this pull request as draft October 28, 2024 18:34
Copy link

@JohnGriffiths JohnGriffiths left a comment

Choose a reason for hiding this comment

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

confirming that the example runs here, and also that it seems to load other data files properly. the latter point I will can make and share a separate NB demonstrating that for 'peace of mind', in lieu of an analysis example discussed on the issue thread. Also noting that given @julien-dubois-k 's comment re some .snirf attributes not being read, so those steps may also need to be added either to this PR or (probably better) a subsequent one. I will be looking at this next.

Afraid I don't know the answer to the uM units Q off hand.

@JohnGriffiths
Copy link

JohnGriffiths commented Nov 7, 2024

@larsoner thanks for this. Have gone through, left review above.

Not quite sure if I missed something specific other than the units Qs and the time_delay_widths Qs that you needed an immediate answer to; LMK

I don't not know to these off the top of my head, but I imagine @julien-dubois-k does. I can take a look at relevant documentation if that is useful.

Re: datasets

LMK if I have got this right re next steps

  • You have two small .snirf files in the example data that is loaded here so new ones of those aren't needed from me now for this test function
  • But you are suggesting to add some other, larger files to mne-misc-data for an analysis example (?)
  • So, would you like me to do a PR to that repo with a .snirf data file for a more fulsome analysis example?
  • If so yes I can get on that. What I would do for the example code is use the colab notebook example @julien-dubois-k linked to in his above comment, and the example file the one used in that nb.
  • ( I can add other finger tapping .snirf files to this if warranted )
  • As I work on that I will take a look at the comments from @julien-dubois-k noted in my review comment above about some .snirf attributes not being read from the .hdf properly by the mne reader. Which may lead to another .io.snirf PR before the analysis example PR.
  • But I think it's prob best not to muddy the current PR with that Q (unless you disagree?)

@julien-dubois-k
Copy link

julien-dubois-k commented Nov 7, 2024 via email

@larsoner
Copy link
Member Author

larsoner commented Nov 7, 2024

You have two small .snirf files in the example data that is loaded here so new ones of those aren't needed from me now for this test function

Well related to the time_delay_widths issue, it would be good to have some small SNIRF test files that actually have this field. I'm not sure if we're using it properly, and if we had test files with this field then we could probably figure it out. @julien-dubois-k if you have such a file we could add it to mne-testing-data!

But you are suggesting to add some other, larger files to mne-misc-data for an analysis example (?) ... What I would do for the example code is use the colab notebook example @julien-dubois-k linked to in his above comment, and the example file the one used in that nb.

Sure, then we modify an example here or more likely in mne-nirs to use it. Doesn't have to be a full-blown example with regression and stats. But if that notebook runs in less than, say, a minute, we could convert it, add some narrative docs if needed, and be good to go. (Examples are meant to run fairly quickly and not require a ton of data to be downloaded!)

( I can add other finger tapping .snirf files to this if warranted )

Sure if it adds something complementary from a documentation perspective. But if it's just additional subjects or a different paradigm from what we currently have, might not be worth it.

(Just an aside, not sure if you've looked at mne-bids-pipeline but there is mne-tools/mne-bids-pipeline#365 that we could resurrect if you're interested in a BIDS pipeline for your data.)

But I think it's prob best not to muddy the current PR with that Q (unless you disagree?)

I don't mind adding more stuff here, a +305/-85 PR is fairly small actually. So based on discussion I think the next steps would be:

  1. @JohnGriffiths open PR to add example data in mne-misc-data that we need for the notebook
  2. I modify mne.datasets code in this PR to reflect the new data
  3. @JohnGriffiths open a PR to mne-nirs to use the new example data based on the notebook
  4. I modify the mne-nirs PR to use this branch temporarily (or you can if you already know how @JohnGriffiths)
  5. @julien-dubois-k supply some tiny files we can add to mne-testing-data that would shed light on the time_delay_widths issue
  6. @julien-dubois-k comment on the unit issue above
  7. I push commits to update tests with new data and/or unit expectations
  8. I make sure all plotting scales etc. are correct based on the MNE-NIRS example
  9. We merge this PR
  10. I remove "use larsoner:feature/TD-nirs_snirf branch" shim from MNE-NIRS PR
  11. We merge the MNE-NIRS PR

@JohnGriffiths
Copy link

@JohnGriffiths open a PR to mne-nirs to use the new example data based on the notebook

Ok I'm on it.

@JohnGriffiths open PR to add example data in mne-misc-data that we need for the notebook

Will do. What is the heuristic on acceptable size for a file in mne-misc-data to use in an example?

I modify the mne-nirs PR to use this branch temporarily (or you can if you already know how @JohnGriffiths)

So just temporarily change this requirements line to the commit url for the current PR, right?

(Just an aside, not sure if you've looked at mne-bids-pipeline but there is mne-tools/mne-bids-pipeline#365 that we could resurrect if you're interested in a BIDS pipeline for your data.)

That's great to know. thanks. We are actually currently in the process of getting an mne-bids-pipeline auto QC in XNAT for our EEG data. Would like to do the same for fNIRS. That would slot in nicely.

@larsoner
Copy link
Member Author

larsoner commented Nov 7, 2024

Will do. What is the heuristic on acceptable size for a file in mne-misc-data to use in an example?

We try to keep our datasets to < 1GB if possible. Looks like mne-misc-data is around 700 MB so it would be good not to add too much more. Maybe up to 100MB is okay for the purposes here? If you need more than that then we can add it as a new dataset on OSF or somewhere else, and add a mne.datasets.<something>.data_path() downloader. That's a little more work but can be worth it.

So just temporarily change this requirements line to the commit url for the current PR, right?

CircleCI is the CI that builds the doc so I would do it here actually

https://github.com/mne-tools/mne-nirs/blob/5c63e855b60d6f7aab907477e582acb323a69d81/tools/circleci_dependencies.sh#L13

@@ -60,6 +62,8 @@
fnirs_fd_ac_amplitude="V",
fnirs_fd_phase="rad",
fnirs_od="V",
fnirs_td_gated_amplitude="M",
Copy link

@julien-dubois-k julien-dubois-k Nov 12, 2024

Choose a reason for hiding this comment

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

these are photon counts, should they be "AU"?

@@ -117,6 +123,8 @@
fnirs_fd_ac_amplitude=1.0,
fnirs_fd_phase=1.0,
fnirs_od=1.0,
fnirs_td_gated_amplitude=1e6,

Choose a reason for hiding this comment

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

I imagine this should be 1.0; same for fnirs_td_moments_amplitude

@@ -60,6 +62,8 @@
fnirs_fd_ac_amplitude="V",
fnirs_fd_phase="rad",
fnirs_od="V",
fnirs_td_gated_amplitude="M",
fnirs_td_moments_amplitude="M",

Choose a reason for hiding this comment

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

this is a bit weird, the dataType Time domain – Moments (TD Moments) - Amplitude can be any of the moments, which all have different units:

  • sum or intensity is photon counts
  • mean is s
  • variance is s^2

Copy link
Member Author

Choose a reason for hiding this comment

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

Ahhh... in this case we probably need to update the FIF "coil type" constants first since these are really sensing/detecting different sorts of data. Like we'll probably need:

  • td_moments_intensity (photon counts)
  • td_moments_mean (s)
  • td_moments_variance (s^2)

or similar

@@ -88,6 +92,8 @@
fnirs_fd_ac_amplitude="V",
fnirs_fd_phase="rad",
fnirs_od="V",
fnirs_td_gated_amplitude="µM",

Choose a reason for hiding this comment

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

update

@@ -57,7 +57,9 @@
)

# Kernel
kernel_hb = testing_path / "SNIRF" / "Kernel" / "Flow50" / "Portal_2021_11" / "hb.snirf"
kernel_path = testing_path / "SNIRF" / "Kernel" / "Flow50" / "Portal_2021_11"

Choose a reason for hiding this comment

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

Here are some files to add. The path could be updated if it matters, these are no longer Flow50/Flow1 but Flow2, and the version is latest (aka, Portal_2024-10-23).

  • Type 201 (Time Domain Gated): c345d04_2.snirf (there is no current equivalent)
  • Type 301 (Time Domain Moments): c345d04_3.snirf (replace td_moments.snirf)
  • Type 99999 (Time Domain Hb Moments): c345d04_5.snirf (replace hb.snirf)

@larsoner
Copy link
Member Author

I am a bit busy for the next two weeks but then should be able to come back to this!

@larsoner larsoner modified the milestones: 1.9, 1.10 Dec 9, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
Development

Successfully merging this pull request may close these issues.

5 participants