From f4a2390b8dcd0281d23bfeb555b2561f6cb6cede Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 30 May 2024 09:59:48 -0400 Subject: [PATCH 1/4] RF: Simplify high-pass filtering in algorithms.confounds Legendre and cosine detrending are implemented almost identically, although with several minor variations. Here I separate regressor creation from detrending to unify the implementations. This now uses `np.linalg.pinv(X)` to estimate the betas in both cases, rather than using `np.linalg.lstsq` in the cosine filter. lstsq uses SVD and can thus fail to converge in rare cases. Under no circumstances should (X.T @ X) be singular, so the pseudoinverse is unique and precisely what we want. --- nipype/algorithms/confounds.py | 84 +++++++++++++++++++--------------- 1 file changed, 46 insertions(+), 38 deletions(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index 4a74ff5cec..4d82bcf080 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -1185,31 +1185,61 @@ def is_outlier(points, thresh=3.5): return timepoints_to_discard -def cosine_filter( - data, timestep, period_cut, remove_mean=True, axis=-1, failure_mode="error" +def _make_cosine_regressors(ntimepoints, timestep, period_cut): + return _full_rank(_cosine_drift(period_cut, timestep * np.arange(ntimepoints)))[0] + + +def _make_legendre_regressors(ntimepoints, degree): + X = np.ones((ntimepoints, 1)) # quick way to calc degree 0 + for i in range(degree): + polynomial_func = Legendre.basis(i + 1) + value_array = np.linspace(-1, 1, ntimepoints) + X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis])) + return X + + +def _high_pass_filter( + data, *args, filter_type, remove_mean=True, axis=-1, failure_mode="error" ): datashape = data.shape timepoints = datashape[axis] if datashape[0] == 0 and failure_mode != "error": return data, np.array([]) - data = data.reshape((-1, timepoints)) + filters = { + 'cosine': _make_cosine_regressors, + 'legendre': _make_legendre_regressors, + } - frametimes = timestep * np.arange(timepoints) - X = _full_rank(_cosine_drift(period_cut, frametimes))[0] + X = filters[filter_type](timepoints, *args) non_constant_regressors = X[:, :-1] if X.shape[1] > 1 else np.array([]) - betas = np.linalg.lstsq(X, data.T)[0] + data = data.reshape((-1, timepoints)) + + betas = np.linalg.pinv(X) @ data.T if not remove_mean: X = X[:, :-1] betas = betas[:-1] - residuals = data - X.dot(betas).T - + residuals = data - (X @ betas).T return residuals.reshape(datashape), non_constant_regressors +def cosine_filter( + data, timestep, period_cut, remove_mean=True, axis=-1, failure_mode="error" +): + return _high_pass_filter( + data, + timestep, + period_cut, + filter_type='cosine', + remove_mean=remove_mean, + axis=axis, + failure_mode=failure_mode, + ) + + def regress_poly(degree, data, remove_mean=True, axis=-1, failure_mode="error"): """ Returns data with degree polynomial regressed out. @@ -1221,36 +1251,14 @@ def regress_poly(degree, data, remove_mean=True, axis=-1, failure_mode="error"): IFLOGGER.debug( "Performing polynomial regression on data of shape %s", str(data.shape) ) - - datashape = data.shape - timepoints = datashape[axis] - if datashape[0] == 0 and failure_mode != "error": - return data, np.array([]) - - # Rearrange all voxel-wise time-series in rows - data = data.reshape((-1, timepoints)) - - # Generate design matrix - X = np.ones((timepoints, 1)) # quick way to calc degree 0 - for i in range(degree): - polynomial_func = Legendre.basis(i + 1) - value_array = np.linspace(-1, 1, timepoints) - X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis])) - - non_constant_regressors = X[:, :-1] if X.shape[1] > 1 else np.array([]) - - # Calculate coefficients - betas = np.linalg.pinv(X).dot(data.T) - - # Estimation - if remove_mean: - datahat = X.dot(betas).T - else: # disregard the first layer of X, which is degree 0 - datahat = X[:, 1:].dot(betas[1:, ...]).T - regressed_data = data - datahat - - # Back to original shape - return regressed_data.reshape(datashape), non_constant_regressors + return _high_pass_filter( + data, + degree, + filter_type='legendre', + remove_mean=remove_mean, + axis=axis, + failure_mode=failure_mode, + ) def combine_mask_files(mask_files, mask_method=None, mask_index=None): From 4717e23b72c1fdef45cb02fa5401bfe7ff45e662 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 30 May 2024 10:46:02 -0400 Subject: [PATCH 2/4] RF: Avoid repeated linspace and hstack calls --- nipype/algorithms/confounds.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index 4d82bcf080..a93d36aca0 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -1190,12 +1190,13 @@ def _make_cosine_regressors(ntimepoints, timestep, period_cut): def _make_legendre_regressors(ntimepoints, degree): - X = np.ones((ntimepoints, 1)) # quick way to calc degree 0 - for i in range(degree): - polynomial_func = Legendre.basis(i + 1) - value_array = np.linspace(-1, 1, ntimepoints) - X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis])) - return X + support = np.linspace([-1], [1], ntimepoints) + return np.hstack( + [ + np.ones((ntimepoints, 1)), # Degree 0 is a constant + *(Legendre.basis(i + 1)(support) for i in range(degree)), + ] + ) def _high_pass_filter( From 510771458fb4e3ad6b3f5c13032c3adc87fe4870 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 30 May 2024 11:11:04 -0400 Subject: [PATCH 3/4] FIX: Normalize constant column to 0th For Legendre regressors, the ith column is the ith-order polynomial, so the constant column is 0. For the cosine regressors, a constant column was appended to the end, in contradiction of the docstring. This brings both into alignment so columns are sorted from lowest to highest frequency and aligns the DCT behavior with its docstring. --- nipype/algorithms/confounds.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index a93d36aca0..ccb22b5066 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -1220,8 +1220,8 @@ def _high_pass_filter( betas = np.linalg.pinv(X) @ data.T if not remove_mean: - X = X[:, :-1] - betas = betas[:-1] + X = X[:, 1:] + betas = betas[1:] residuals = data - (X @ betas).T return residuals.reshape(datashape), non_constant_regressors @@ -1547,9 +1547,9 @@ def _cosine_drift(period_cut, frametimes): nfct = np.sqrt(2.0 / len_tim) for k in range(1, order): - cdrift[:, k - 1] = nfct * np.cos((np.pi / len_tim) * (n_times + 0.5) * k) + cdrift[:, k] = nfct * np.cos((np.pi / len_tim) * (n_times + 0.5) * k) - cdrift[:, order - 1] = 1.0 # or 1./sqrt(len_tim) to normalize + cdrift[:, 0] = 1.0 # or 1./sqrt(len_tim) to normalize return cdrift From 17bac0817a36d9430ec25e3f68ec80379ffa9752 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 30 May 2024 11:21:15 -0400 Subject: [PATCH 4/4] RF: Slightly simplify cosine implementation --- nipype/algorithms/confounds.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index ccb22b5066..dba0545e8c 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -1535,7 +1535,6 @@ def _cosine_drift(period_cut, frametimes): Ref: http://en.wikipedia.org/wiki/Discrete_cosine_transform DCT-II """ len_tim = len(frametimes) - n_times = np.arange(len_tim) hfcut = 1.0 / period_cut # input parameter is the period # frametimes.max() should be (len_tim-1)*dt @@ -1543,13 +1542,15 @@ def _cosine_drift(period_cut, frametimes): # hfcut = 1/(2*dt) yields len_time # If series is too short, return constant regressor order = max(int(np.floor(2 * len_tim * hfcut * dt)), 1) - cdrift = np.zeros((len_tim, order)) - nfct = np.sqrt(2.0 / len_tim) + cdrift = np.ones((len_tim, order)) - for k in range(1, order): - cdrift[:, k] = nfct * np.cos((np.pi / len_tim) * (n_times + 0.5) * k) + if order > 1: + nfct = np.sqrt(2.0 / len_tim) + support = (np.pi / len_tim) * (np.arange(len_tim) + 0.5) + + for k in range(1, order): + cdrift[:, k] = nfct * np.cos(support * k) - cdrift[:, 0] = 1.0 # or 1./sqrt(len_tim) to normalize return cdrift