Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Sign up
Appearance settings

Commit f4a2390

Browse files
committed
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.
1 parent 4d1352a commit f4a2390

File tree

1 file changed

+46
-38
lines changed

1 file changed

+46
-38
lines changed

‎nipype/algorithms/confounds.py‎

Lines changed: 46 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -1185,31 +1185,61 @@ def is_outlier(points, thresh=3.5):
11851185
return timepoints_to_discard
11861186

11871187

1188-
def cosine_filter(
1189-
data, timestep, period_cut, remove_mean=True, axis=-1, failure_mode="error"
1188+
def _make_cosine_regressors(ntimepoints, timestep, period_cut):
1189+
return _full_rank(_cosine_drift(period_cut, timestep * np.arange(ntimepoints)))[0]
1190+
1191+
1192+
def _make_legendre_regressors(ntimepoints, degree):
1193+
X = np.ones((ntimepoints, 1)) # quick way to calc degree 0
1194+
for i in range(degree):
1195+
polynomial_func = Legendre.basis(i + 1)
1196+
value_array = np.linspace(-1, 1, ntimepoints)
1197+
X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis]))
1198+
return X
1199+
1200+
1201+
def _high_pass_filter(
1202+
data, *args, filter_type, remove_mean=True, axis=-1, failure_mode="error"
11901203
):
11911204
datashape = data.shape
11921205
timepoints = datashape[axis]
11931206
if datashape[0] == 0 and failure_mode != "error":
11941207
return data, np.array([])
11951208

1196-
data = data.reshape((-1, timepoints))
1209+
filters = {
1210+
'cosine': _make_cosine_regressors,
1211+
'legendre': _make_legendre_regressors,
1212+
}
11971213

1198-
frametimes = timestep * np.arange(timepoints)
1199-
X = _full_rank(_cosine_drift(period_cut, frametimes))[0]
1214+
X = filters[filter_type](timepoints, *args)
12001215
non_constant_regressors = X[:, :-1] if X.shape[1] > 1 else np.array([])
12011216

1202-
betas = np.linalg.lstsq(X, data.T)[0]
1217+
data = data.reshape((-1, timepoints))
1218+
1219+
betas = np.linalg.pinv(X) @ data.T
12031220

12041221
if not remove_mean:
12051222
X = X[:, :-1]
12061223
betas = betas[:-1]
12071224

1208-
residuals = data - X.dot(betas).T
1209-
1225+
residuals = data - (X @ betas).T
12101226
return residuals.reshape(datashape), non_constant_regressors
12111227

12121228

1229+
def cosine_filter(
1230+
data, timestep, period_cut, remove_mean=True, axis=-1, failure_mode="error"
1231+
):
1232+
return _high_pass_filter(
1233+
data,
1234+
timestep,
1235+
period_cut,
1236+
filter_type='cosine',
1237+
remove_mean=remove_mean,
1238+
axis=axis,
1239+
failure_mode=failure_mode,
1240+
)
1241+
1242+
12131243
def regress_poly(degree, data, remove_mean=True, axis=-1, failure_mode="error"):
12141244
"""
12151245
Returns data with degree polynomial regressed out.
@@ -1221,36 +1251,14 @@ def regress_poly(degree, data, remove_mean=True, axis=-1, failure_mode="error"):
12211251
IFLOGGER.debug(
12221252
"Performing polynomial regression on data of shape %s", str(data.shape)
12231253
)
1224-
1225-
datashape = data.shape
1226-
timepoints = datashape[axis]
1227-
if datashape[0] == 0 and failure_mode != "error":
1228-
return data, np.array([])
1229-
1230-
# Rearrange all voxel-wise time-series in rows
1231-
data = data.reshape((-1, timepoints))
1232-
1233-
# Generate design matrix
1234-
X = np.ones((timepoints, 1)) # quick way to calc degree 0
1235-
for i in range(degree):
1236-
polynomial_func = Legendre.basis(i + 1)
1237-
value_array = np.linspace(-1, 1, timepoints)
1238-
X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis]))
1239-
1240-
non_constant_regressors = X[:, :-1] if X.shape[1] > 1 else np.array([])
1241-
1242-
# Calculate coefficients
1243-
betas = np.linalg.pinv(X).dot(data.T)
1244-
1245-
# Estimation
1246-
if remove_mean:
1247-
datahat = X.dot(betas).T
1248-
else: # disregard the first layer of X, which is degree 0
1249-
datahat = X[:, 1:].dot(betas[1:, ...]).T
1250-
regressed_data = data - datahat
1251-
1252-
# Back to original shape
1253-
return regressed_data.reshape(datashape), non_constant_regressors
1254+
return _high_pass_filter(
1255+
data,
1256+
degree,
1257+
filter_type='legendre',
1258+
remove_mean=remove_mean,
1259+
axis=axis,
1260+
failure_mode=failure_mode,
1261+
)
12541262

12551263

12561264
def combine_mask_files(mask_files, mask_method=None, mask_index=None):

0 commit comments

Comments
(0)

AltStyle によって変換されたページ (->オリジナル) /