From 4c49e9061e3c7c36e187cccddc61ff4d33e4ead2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eneko=20Uru=C3=B1uela?= Date: Sat, 20 Dec 2025 10:26:27 -0700 Subject: [PATCH 1/3] Update 05_3dMEPFM.md for pySPFM v2.0 API MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Migrate notebook from legacy pySPFM.pySPFM() function to new scikit-learn style SparseDeconvolution class API: - Use SparseDeconvolution with fit/transform pattern - Add NiftiMasker for data loading and inverse transformation - Concatenate multi-echo data along time axis as required by new API - Add documentation of model attributes (coef_, lambda_, hrf_matrix_) - Save fitted signal and residuals as additional outputs 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 --- content/05_3dMEPFM.md | 65 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 58 insertions(+), 7 deletions(-) diff --git a/content/05_3dMEPFM.md b/content/05_3dMEPFM.md index 6ed815e..eb97223 100644 --- a/content/05_3dMEPFM.md +++ b/content/05_3dMEPFM.md @@ -17,7 +17,9 @@ import json import os from glob import glob -import nibabel as nb +import nibabel as nib +import numpy as np +from nilearn.maskers import NiftiMasker data_path = os.path.abspath('../DATA') ``` @@ -53,16 +55,65 @@ out_dir = os.path.join(data_path, "pySPFM") ```{code-cell} ipython3 :tags: [output_scroll] -from pySPFM import pySPFM +from pySPFM import SparseDeconvolution -pySPFM.pySPFM( - data_fn=data_files, - mask_fn=mask_file, - output_filename=os.path.join(out_dir, "out"), +# Create masker to convert 4D NIfTI data to 2D array +masker = NiftiMasker(mask_img=mask_file) + +# Load and mask each echo, then concatenate along time axis +# For multi-echo data, timepoints from different echoes are concatenated along the first axis +masked_data = [] +for f in data_files: + echo_data = masker.fit_transform(f) # Shape: (n_timepoints, n_voxels) + masked_data.append(echo_data) + +X = np.vstack(masked_data) # Shape: (n_timepoints * n_echoes, n_voxels) + +# Fit the sparse deconvolution model +model = SparseDeconvolution( tr=2.47, - out_dir=out_dir, te=echo_times, + criterion="bic", ) +model.fit(X) + +# Get the deconvolved activity-inducing signals +activity = model.coef_ # Shape: (n_timepoints, n_voxels) + +# Transform back to NIfTI image and save +os.makedirs(out_dir, exist_ok=True) +activity_img = masker.inverse_transform(activity) +activity_img.to_filename(os.path.join(out_dir, "out_activity.nii.gz")) + +# Also save the regularization parameter values +np.save(os.path.join(out_dir, "out_lambda.npy"), model.lambda_) + +print(f"Activity shape: {activity.shape}") +print(f"Saved activity to: {os.path.join(out_dir, 'out_activity.nii.gz')}") +``` + +The `SparseDeconvolution` model provides several useful attributes and methods after fitting: + +- `coef_`: The deconvolved activity-inducing signals +- `lambda_`: The regularization parameter values +- `hrf_matrix_`: The HRF convolution matrix used +- `get_fitted_signal()`: Returns the fitted (reconstructed) signal +- `get_residuals(X)`: Returns the residuals between the original data and fitted signal + +```{code-cell} ipython3 +# Get the fitted signal and residuals +fitted_signal = model.get_fitted_signal() +residuals = model.get_residuals(X) + +# Save additional outputs +fitted_img = masker.inverse_transform(fitted_signal) +fitted_img.to_filename(os.path.join(out_dir, "out_fitted.nii.gz")) + +residuals_img = masker.inverse_transform(residuals) +residuals_img.to_filename(os.path.join(out_dir, "out_residuals.nii.gz")) + +print(f"Fitted signal shape: {fitted_signal.shape}") +print(f"Residuals shape: {residuals.shape}") ``` The pySPFM workflow writes out a number of files. From d8719b35882b3184f09e38a45f0aefc3f1049e9b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eneko=20Uru=C3=B1uela?= Date: Sat, 20 Dec 2025 11:00:19 -0700 Subject: [PATCH 2/3] Address PR review comments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Remove unused nibabel import - Fit NiftiMasker once before loop, use transform() inside - Clarify comments about multi-echo data stacking - Add shape information to coef_ documentation - Clarify method parameter requirements in documentation - Save fitted_signal and residuals as numpy arrays instead of NIfTI (their shape doesn't match single-echo masker expectations) 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 --- content/05_3dMEPFM.md | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/content/05_3dMEPFM.md b/content/05_3dMEPFM.md index eb97223..d98f6b6 100644 --- a/content/05_3dMEPFM.md +++ b/content/05_3dMEPFM.md @@ -17,7 +17,6 @@ import json import os from glob import glob -import nibabel as nib import numpy as np from nilearn.maskers import NiftiMasker @@ -60,14 +59,17 @@ from pySPFM import SparseDeconvolution # Create masker to convert 4D NIfTI data to 2D array masker = NiftiMasker(mask_img=mask_file) -# Load and mask each echo, then concatenate along time axis -# For multi-echo data, timepoints from different echoes are concatenated along the first axis +# Fit masker once on a representative image (first echo) +masker.fit(data_files[0]) + +# Load and mask each echo, then concatenate along the time axis +# For multi-echo data, each echo's data (shape: n_timepoints × n_voxels) are stacked sequentially masked_data = [] for f in data_files: - echo_data = masker.fit_transform(f) # Shape: (n_timepoints, n_voxels) + echo_data = masker.transform(f) # Shape: (n_timepoints, n_voxels) masked_data.append(echo_data) -X = np.vstack(masked_data) # Shape: (n_timepoints * n_echoes, n_voxels) +X = np.vstack(masked_data) # Shape: (n_echoes * n_timepoints, n_voxels) # Fit the sparse deconvolution model model = SparseDeconvolution( @@ -78,7 +80,9 @@ model = SparseDeconvolution( model.fit(X) # Get the deconvolved activity-inducing signals -activity = model.coef_ # Shape: (n_timepoints, n_voxels) +# Note: coef_ has shape (n_timepoints, n_voxels) - the model recovers +# the underlying neural activity at the original temporal resolution +activity = model.coef_ # Transform back to NIfTI image and save os.makedirs(out_dir, exist_ok=True) @@ -94,23 +98,22 @@ print(f"Saved activity to: {os.path.join(out_dir, 'out_activity.nii.gz')}") The `SparseDeconvolution` model provides several useful attributes and methods after fitting: -- `coef_`: The deconvolved activity-inducing signals +- `coef_`: The deconvolved activity-inducing signals (shape: n_timepoints × n_voxels) - `lambda_`: The regularization parameter values - `hrf_matrix_`: The HRF convolution matrix used -- `get_fitted_signal()`: Returns the fitted (reconstructed) signal -- `get_residuals(X)`: Returns the residuals between the original data and fitted signal +- `get_fitted_signal()`: Returns the fitted (reconstructed) signal; takes no arguments +- `get_residuals(X)`: Returns the residuals between the original data and fitted signal; requires the input data `X` as an argument ```{code-cell} ipython3 # Get the fitted signal and residuals +# Note: These have shape (n_echoes * n_timepoints, n_voxels) matching the input X fitted_signal = model.get_fitted_signal() residuals = model.get_residuals(X) -# Save additional outputs -fitted_img = masker.inverse_transform(fitted_signal) -fitted_img.to_filename(os.path.join(out_dir, "out_fitted.nii.gz")) - -residuals_img = masker.inverse_transform(residuals) -residuals_img.to_filename(os.path.join(out_dir, "out_residuals.nii.gz")) +# Save the fitted signal and residuals as numpy arrays +# (shape doesn't match single-echo masker expectations for NIfTI output) +np.save(os.path.join(out_dir, "out_fitted.npy"), fitted_signal) +np.save(os.path.join(out_dir, "out_residuals.npy"), residuals) print(f"Fitted signal shape: {fitted_signal.shape}") print(f"Residuals shape: {residuals.shape}") From d359e4992e36c266104574522ed3dd5f02ce4718 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eneko=20Uru=C3=B1uela?= Date: Sat, 20 Dec 2025 14:01:29 -0700 Subject: [PATCH 3/3] Enable parallel processing in SparseDeconvolution MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add n_jobs=-1 to use all available CPU cores, significantly speeding up the voxel-wise deconvolution. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 --- content/05_3dMEPFM.md | 1 + 1 file changed, 1 insertion(+) diff --git a/content/05_3dMEPFM.md b/content/05_3dMEPFM.md index d98f6b6..4db2505 100644 --- a/content/05_3dMEPFM.md +++ b/content/05_3dMEPFM.md @@ -76,6 +76,7 @@ model = SparseDeconvolution( tr=2.47, te=echo_times, criterion="bic", + n_jobs=-1, # Use all available CPU cores ) model.fit(X)