ENH: python: Add magnitude_only mode, fix three bugs, tighten output dtype and docstrings

[MoniDoerig]

Note: issues appear disabled on this repository, so opening this PR directly.

This PR adds support for ARG.magnitude_only=1 mode in the Python NORDIC implementation (python/src/nordic/denoise.py), fixes three bugs in the same file, tightens the output dtype and docstrings, and adds pytest tests covering the magnitude_only mode and dtype behaviour (python/tests/test_matlab_parity.py). It includes an output-comparison notebook on a public dataset (OpenNeuro ds004928) showing the resulting Python output matches MATLAB NIFTI_NORDIC across all three algorithm flavours.

Changes

Bug fixes in python/src/nordic/denoise.py

Three issues in subfunction_loop_for_nvr_avg_update:

1. rangeData length mismatch — vals[: R - centering + 1] is one element longer than rangeMP, breaking the element-wise division rangeData / rangeMP. Fixed by trimming the trailing element.
2. Empty np.where crash — when sigmasq_2 < sigmasq_1 is empty (noise never exceeds signal), the code raised IndexError. Now interpreted as "no components are noise" — keep all singular values.
3. Out-of-bounds noise indexing — companion to (2): first_removed_component = S.size can exceed len(sigmasq_2); guarded in both update paths.

New: magnitude_only mode

Mirrors MATLAB ARG.magnitude_only=1:

• Explicit — magnitude_only=True forces temporal_phase=0 and ignores any pha_file (with warnings).
• Implicit safety net — pha_file=None with temporal_phase != 0 is auto-corrected with a warning. Previously this silently applied a phase filter to non-existent phase data, introducing a spurious complex modulation.
• Behaviour covered by python/tests/test_matlab_parity.py::test_magnitude_only_ (3 tests, synthetic phantom, runs in seconds).

Output dtype

• Force float32 on the magnitude and phase outputs (matches MATLAB info.Datatype = single). Without this, integer-typed inputs (e.g. int16 BOLD) were silently truncated by nibabel on write.
• Behaviour covered by python/tests/test_matlab_parity.py::test_output_is_float32_for_int16_input.

Docstring clarifications

Document the actual defaults for kernel_size_gfactor ([14, 14, 1, 90], with the temporal axis capped at n_vols) and kernel_size_pca (cubic patch with 11:1 spatial:temporal ratio, with slice-clamping when n_slices < cubic edge).

Validation

Tested against the MATLAB reference (NIFTI_NORDIC.m) on three subjects of the publicly available OpenNeuro ds004928 (Faes/Vizioli/Moeller, 7T submillimetre auditory fMRI with magnitude, phase, and trailing noise volumes), across all three documented algorithm flavours.

Headline numbers:

Subject	algorithm	Pearson r	rel-diff median	MATLAB/Python tSNR ratio
sub-01	nordic	0.999974	0.690 %	1.0071
sub-01	gfactor+mppca	0.999999	0.020 %	1.0000
sub-01	mppca	1.000000	0.017 %	1.0000
sub-02	nordic	0.999990	0.134 %	1.0007
sub-02	gfactor+mppca	0.999993	0.010 %	1.0000
sub-02	mppca	1.000000	0.008 %	1.0000
sub-03	nordic	0.999987	0.258 %	1.0025
sub-03	gfactor+mppca	0.999994	0.009 %	0.9999
sub-03	mppca	1.000000	0.007 %	1.0000

All 9 cells clear r ≥ 0.9999 voxel-wise, median rel-diff < 1 %, and tSNR ratio within ~1 % of unity. Residuals scale with algorithmic complexity (more LAPACK-ordering exposure → looser parity), consistent with float32 SVD ordering noise.

The notebook (python/validation/matlab_parity_ds004928.ipynb) contains the full reproducibility recipe — DataLad install through both Python and MATLAB invocations — so the result can be replicated end-to-end from public data.

[tsalo]

Thank you for addressing those bugs- they all look solid. I'm not really a fan of adding a magnitude_only parameter though. I think it's fine to infer if phase data should be used by if phase data are present, even if that doesn't align with the original MATLAB implementation.

[MoniDoerig]

Thanks for the review! The four suggestions are all clear improvements, I've just applied them.

I'm not really a fan of adding a magnitude_only parameter though.

On magnitude_only: My one concern with pure inference is silent fallback — if a caller meant to pass phase data but forgot, they'd silently get a different algorithm with no warning, and Vizioli 2021 shows complex NORDIC removes meaningfully more thermal noise.

Would a compromise work: drop the parameter as you suggest, but emit a UserWarning when pha_file is None and temporal_phase != 0 ("no phase data provided; running magnitude-only NORDIC")? That keeps the inference-from-arguments API and prevents the silent-fallback case. Happy to push that as a follow-up commit.