[FIX] Disable extrapolation of out-of-bounds volumes

[jhuguetn]

• Closes [BUG] Interpolation of high motion volumes at the beginning/end scan #4024

Changes proposed in this pull request:

• include a extrapolate parameter in the signal._interpolate_volumes() function, exposing the identically named parameter from the Scipy class used underneath for the Splines-based interpolation (see https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicSpline.html#scipy.interpolate.CubicSpline).
• setting the default extrapolate new parameter value to False, disabling the undesired interpolation on out-of-bounds volumes

[github-actions]

👋 @jhuguetn Thanks for creating a PR!

Until this PR is ready for review, you can include the [WIP] tag in its title, or leave it as a github draft.

Please make sure it is compliant with our contributing guidelines. In particular, be sure it checks the boxes listed below.

• PR has an interpretable title.
• PR links to Github issue with mention Closes #XXXX (see our documentation on PR structure)
• Code is PEP8-compliant (see our documentation on coding style)
• Changelog or what's new entry in doc/changes/latest.rst (see our documentation on PR structure)

For new features:

• There is at least one unit test per new function / class (see our documentation on testing)
• The new feature is demoed in at least one relevant example.

For bug fixes:

• There is at least one test that would fail under the original bug conditions.

We will review it as quick as possible, feel free to ping us with questions if needed.

[ymzayek]

Hi @jhuguetn thanks for opening the PR! Could you also add a test to nilearn/tests/test_signal.py?
I'm also wondering if this needs to exposed in the public function signal.clean that ends up calling this function? If there is a use case for that then it would be possible to pass it as a butterworth kwarg but then it needs to be added as a parameter to _handle_scrubbed_volumes.

[codecov]

Codecov Report

Merging #4028 (bd8e627) into main (912e239) will decrease coverage by 0.08%.
Report is 40 commits behind head on main.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main    #4028      +/-   ##
==========================================
- Coverage   91.59%   91.52%   -0.08%     
==========================================
  Files         145      143       -2     
  Lines       16238    16086     -152     
  Branches     3376     3342      -34     
==========================================
- Hits        14874    14722     -152     
- Misses        817      819       +2     
+ Partials      547      545       -2     

Flag	Coverage Δ
macos-latest_3.10	91.51% <100.00%> (-0.01%)	⬇️
macos-latest_3.11	91.51% <100.00%> (-0.01%)	⬇️
macos-latest_3.12	91.51% <100.00%> (-0.01%)	⬇️
macos-latest_3.8	?
macos-latest_3.9	91.48% <100.00%> (-0.01%)	⬇️
ubuntu-latest_3.10	91.51% <100.00%> (-0.01%)	⬇️
ubuntu-latest_3.11	91.51% <100.00%> (-0.01%)	⬇️
ubuntu-latest_3.12	91.51% <100.00%> (-0.01%)	⬇️
ubuntu-latest_3.8	?
ubuntu-latest_3.9	91.48% <100.00%> (-0.01%)	⬇️
windows-latest_3.10	91.47% <100.00%> (+<0.01%)	⬆️
windows-latest_3.11	91.47% <100.00%> (+<0.01%)	⬆️
windows-latest_3.12	91.47% <100.00%> (+<0.01%)	⬆️
windows-latest_3.8	?
windows-latest_3.9	?

Flags with carried forward coverage won't be shown. Click here to find out more.

Files	Coverage Δ
nilearn/signal.py	97.21% <100.00%> (+0.06%)	⬆️

... and 77 files with indirect coverage changes

📣 Codecov offers a browser extension for seamless coverage viewing on GitHub. Try it in Chrome or Firefox today!

[htwangtw]

This change should be reflected in the parameter of signal.clean.
I raised a few things and would like @ymzayek and @Remi-Gau to have a say.
Also we need a test for this change

[Remi-Gau]

Could you also add a test to nilearn/tests/test_signal.py?

I agree that we definitely need a test for (preferably not just a smoke test)

[ymzayek]

Initially it seemed that this is a bug fix with limited scope, hence the current solution could be ok. However, it is a behavioral change and to me the impact on user code is hard to tell so I agree with @htwangtw that it is safest to go through a deprecation cycle especially since the faulty implementation did not actually fail any user code IIUC. So I would follow the suggestions as laid out by @htwangtw and @Remi-Gau to expose extrapolate in signal clean, keep default behavior, throw a deprecation warning in the appropriate place to tell the user to set extrapolate=False. The warning should say that the default behavior will change in release 0.13. @jhuguetn is it clear?

[jhuguetn]

Yes, it is @ymzayek. Thanks!

[bthirion]

You probably want to add a test to ensure that no extrapolation is performed.

[jhuguetn]

You probably want to add a test to ensure that no extrapolation is performed.

Yes, I am working on that.

I have encountered a small bump that perhaps any of you all can help me with. When extrapolation is disabled, scipy CubicSpline returns NaNs for those out-of-bounds data points (i.e. volumes to be censored at the ends of the time series). See below example with n=5 first data points being censored:

While testing the new extrapolate switch parameter inside the signal.clean() function -as complex as it is- I've encountered that other steps in the signal.clean() workflow do not expect having NaNs in both the signal/confounds discarded (due to avoiding extrapolation) and raise errors:

# Pivoting in qr decomposition was added in scipy 0.10
Q, R, _ = linalg.qr(confounds, mode="economic", pivoting=True)
...
ValueError: array must not contain infs or NaNs

To cope with that I've tried dropping the NaNs, but then by design some posterior steps in charge of volume censoring in signal.clean() complain as well since censored volumes/samples have already been suppressed beforehand:

def _censor_signals(signals, confounds, sample_mask):
      """Apply sample masks to data."""
      signals = signals[sample_mask, :]
...      
IndexError: index 35 is out of bounds for axis 0 with size 35

As an alternative, I thought we could assign a fix/hardcoded value temporarily to those NaNs (which would be censored later) or replace them by the original data, but still I don't like the idea since I lose notion of what's been discarded in the interpolation step. In addition, it might have an impact doing so on posterior steps. As an example see below a signal.clean() test execution replacing the NaNs by zeroes,

interpolated_signal = clean(signals, confounds=confounds, sample_mask=sample_mask, filter='butterworth', 
                            low_pass=0.08, high_pass=0.01, t_r=2.5, extrapolate=False)

I am interested in knowing what you folks think. Thanks!

[htwangtw]

The NaN issue mentioned by @jhuguetn - I pointed out this potential issue in the original discussion. I could have elaborated it more.

After the _handle_scrubbed_volumes step, you should get signals and confounds the same shape, with the heads and tails trimmed if need be.
You can reuse the _censor_signals function to remove NaN after this block:

nilearn/nilearn/signal.py

Lines 821 to 826 in db5a970

	if filter_type == "butterworth":
	signals = _interpolate_volumes(signals, sample_mask, t_r,
	extrapolate)
	if confounds is not None:
	confounds = _interpolate_volumes(confounds, sample_mask, t_r,
	extrapolate)

What you will need a new sample_mask that's for the volumes of the heads and tails you want to discard.
Let me know if this makes sense. I can implement it if need be.

[jhuguetn]

Thanks @htwangtw for the suggestion, I will try that

[ymzayek]

Hi @jhuguetn did you accidentally close this or would you like one of the core developers to take over working on it?

[jhuguetn]

Sorry, I discarded the first changes and replaced them with new ones based on suggestions made.

[jhuguetn]

I have now included a proposal for dealing with the NaNs issue following your advice @htwangtw of reseting the sample_mask indexes. Thanks!

After the _handle_scrubbed_volumes step, you should get signals and confounds the same shape, with the heads and tails trimmed if need be. You can reuse the _censor_signals function to remove NaN

The problem of using _censor_signals function here is that not only suppresses the volumes not interpolated in the ends of the data as desired but also other ones that might be in the middle (which were interpolated and you actually want include in the filtering step!). Instead, I propose simply removing the NaNs from interpolated volumes and resetting the indexes when volumes at the very beginning are censored.

In addition, my PR includes now a test for _handle_scrubbed_volumes asserting the differences in using or not extrapolation in terms of data included in the interpolation substep. Perhaps we could also include a more general test for clean?

[htwangtw]

We need a new test for clean as there are code modified beyond the private function.

I suggest to look for all instance of sample_mask after _handle_scrubbed_volumes and see what's the impact of the NaN values.

Let me know if you would prefer me to take over.

[bthirion]

We're almost there I think ! Thx,

[bthirion]

LGTM, thx.

[ymzayek]

@jhuguetn can you add a test checking that the warning is raised? And can you also add a whatsnew entry to doc/changes/latest.rst and your name to CITATION.cff? (see https://nilearn.github.io/dev/development.html#changelog for more detail)

[htwangtw]

I would recommand adding one more test in nilearn/tests/test_signal.py::test_sample_mask to test all the subsequent processes are performed correctly.

[jhuguetn]

@jhuguetn can you add a test checking that the warning is raised? And can you also add a whatsnew entry to doc/changes/latest.rst and your name to CITATION.cff? (see https://nilearn.github.io/dev/development.html#changelog for more detail)

Done, see 41635d4 and 8d4e444.

[htwangtw]

LGTM now! Thanks for all the hard work!

[ymzayek]

Thx @jhuguetn ! The failure is unrelated to this PR. Merging