From 2c31a4bf3bd6358c7a9bba468b836b4010bc6163 Mon Sep 17 00:00:00 2001 From: Eduard Ort Date: Thu, 2 Dec 2021 15:22:21 +0100 Subject: [PATCH 1/6] memsaving suggestions --- meegkit/dss.py | 29 ++++++++++++----------------- meegkit/tspca.py | 13 ++++++------- meegkit/utils/denoise.py | 13 +++++++++++-- meegkit/utils/matrix.py | 5 +++-- meegkit/utils/sig.py | 10 +++++----- 5 files changed, 37 insertions(+), 33 deletions(-) diff --git a/meegkit/dss.py b/meegkit/dss.py index b67b81b7..87e915c5 100644 --- a/meegkit/dss.py +++ b/meegkit/dss.py @@ -197,24 +197,21 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, blocksize = n_samples # Recentre data - X = demean(X) + X = demean(X, inplace=True) # Cancel line_frequency and harmonics + light lowpass X_filt = smooth(X, sfreq / fline) - - # Subtract clean data from original data. The result is the artifact plus - # some residual biological signal - X_noise = X - X_filt - + + # X - X_filt results in the artifact plus some residual biological signal # Reduce dimensionality to avoid overfitting if nkeep is not None: - cov_X_res = tscov(X_noise)[0] + cov_X_res = tscov(X - X_filt)[0] V, _ = pca(cov_X_res, nkeep) - X_noise_pca = X_noise @ V + X_noise_pca = (X - X_filt) @ V else: - X_noise_pca = X_noise.copy() + X_noise_pca = (X - X_filt).copy() nkeep = n_chans - + # Compute blockwise covariances of raw and biased data n_harm = np.floor((sfreq / 2) / fline).astype(int) c0 = np.zeros((nkeep, nkeep)) @@ -226,10 +223,9 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, X_block = X_block.transpose(1, 2, 0) # bias data - X_bias = gaussfilt(X_block, sfreq, fline, fwhm=1, n_harm=n_harm) c0 += tscov(X_block)[0] - c1 += tscov(X_bias)[0] - + c1 += tscov(gaussfilt(X_block, sfreq, fline, fwhm=1, n_harm=n_harm))[0] + # DSS to isolate line components from residual todss, _, pwr0, pwr1 = dss0(c0, c1) @@ -244,16 +240,15 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, # Remove line components from X_noise idx_remove = np.arange(nremove) X_artifact = matmul3d(X_noise_pca, todss[:, idx_remove]) - X_res = tsr(X_noise, X_artifact)[0] # project them out - + X_res = tsr(X - X_filt, X_artifact)[0] # project them out # reconstruct clean signal y = X_filt + X_res - artifact = X - y # Power of components p = wpwr(X - y)[0] / wpwr(X)[0] print('Power of components removed by DSS: {:.2f}'.format(p)) - return y, artifact + # return the reconstructed clean signal, and the artifact + return y, X - y def dss_line_iter(data, fline, sfreq, win_sz=10, spot_sz=2.5, diff --git a/meegkit/tspca.py b/meegkit/tspca.py index 9769bffc..81ee64aa 100644 --- a/meegkit/tspca.py +++ b/meegkit/tspca.py @@ -107,6 +107,7 @@ def tsr(X, R, shifts=None, wX=None, wR=None, keep=None, thresh=1e-12): Weights applied by TSR. """ + from psutil import Process ndims = X.ndim X = unsqueeze(X) R = unsqueeze(R) @@ -158,8 +159,8 @@ def tsr(X, R, shifts=None, wX=None, wR=None, keep=None, thresh=1e-12): wR = weights # remove weighted means - X, mean1 = demean(X, wX, return_mean=True) - R = demean(R, wR) + X, mean1 = demean(X, wX, return_mean=True, inplace=True) + R = demean(R, wR, inplace=True) # equalize power of R channels, the equalize power of the R PCs # if R.shape[1] > 1: @@ -182,11 +183,9 @@ def tsr(X, R, shifts=None, wX=None, wR=None, keep=None, thresh=1e-12): y = np.zeros((n_samples_X, n_chans_X, n_trials_X)) for t in np.arange(n_trials_X): r = multishift(R[..., t], shifts, reshape=True) - z = r @ regression - y[..., t] = X[:z.shape[0], :, t] - z - - y, mean2 = demean(y, wX, return_mean=True) - + y[..., t] = X[:z.shape[0], :, t] - (r @ regression) + + y, mean2 = demean(y, wX, return_mean=True, inplace=True) idx = np.arange(offset1, initial_samples - offset2) mean_total = mean1 + mean2 weights = wR diff --git a/meegkit/utils/denoise.py b/meegkit/utils/denoise.py index dc87bcd7..9181397b 100644 --- a/meegkit/utils/denoise.py +++ b/meegkit/utils/denoise.py @@ -7,7 +7,7 @@ from .matrix import fold, theshapeof, unfold, _check_weights -def demean(X, weights=None, return_mean=False): +def demean(X, weights=None, return_mean=False, inplace=False): """Remove weighted mean over rows (samples). Parameters @@ -15,6 +15,9 @@ def demean(X, weights=None, return_mean=False): X : array, shape=(n_samples, n_channels[, n_trials]) Data. weights : array, shape=(n_samples) + return_mean : bool, optional + inplace : bool, optional + Save the resulting array in X, or in a new array Returns ------- @@ -43,9 +46,15 @@ def demean(X, weights=None, return_mean=False): raise ValueError('Weight array should have either the same ' + 'number of columns as X array, or 1 column.') - demeaned_X = X - mn else: mn = np.mean(X, axis=0, keepdims=True) + + # return new array or remove the mean in X itself + if inplace: + np.subtract(X, mn, out=X) + # create alias, so that rest of script can run with the same variable + demeaned_X = X + else: demeaned_X = X - mn if n_trials > 1 or ndims == 3: diff --git a/meegkit/utils/matrix.py b/meegkit/utils/matrix.py index bacc3580..c7aa85fe 100644 --- a/meegkit/utils/matrix.py +++ b/meegkit/utils/matrix.py @@ -503,8 +503,9 @@ def fold(X, epoch_size): n_chans = X.shape[0] // epoch_size if X.shape[0] / epoch_size >= 1: - X = np.transpose(np.reshape(X, (epoch_size, n_chans, X.shape[1]), - order="F").copy(), [0, 2, 1]) + X = np.reshape(X, (epoch_size, X.shape[1], n_chans), order="F") + #X = np.transpose(np.reshape(X, (epoch_size, n_chans, X.shape[1]), + # order="F").copy(), [0, 2, 1]) return X diff --git a/meegkit/utils/sig.py b/meegkit/utils/sig.py index c8a312fe..499ca272 100644 --- a/meegkit/utils/sig.py +++ b/meegkit/utils/sig.py @@ -345,13 +345,13 @@ def gaussfilt(data, srate, f, fwhm, n_harm=1, shift=0, return_empvals=False, fx = fx + gauss # filter - tmp = np.fft.fft(data, axis=0) if data.ndim == 2: - tmp *= fx[:, None] + filtdat = 2 * np.real(np.fft.ifft( + np.fft.fft(data, axis=0)*fx[:, None], axis=0)) elif data.ndim == 3: - tmp *= fx[:, None, None] - - filtdat = 2 * np.real(np.fft.ifft(tmp, axis=0)) + filtdat = 2 * np.real(np.fft.ifft( + np.fft.fft(data, axis=0)*fx[:, None, None], + axis=0)) if return_empvals or show: empVals[0] = hz[idx_p] From 6dbd83c77e3e7b203c46b1aceed1b438ad93f787 Mon Sep 17 00:00:00 2001 From: Eduard Ort Date: Thu, 2 Dec 2021 15:24:58 +0100 Subject: [PATCH 2/6] remove testing code --- meegkit/tspca.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/meegkit/tspca.py b/meegkit/tspca.py index 81ee64aa..0efbac21 100644 --- a/meegkit/tspca.py +++ b/meegkit/tspca.py @@ -107,7 +107,6 @@ def tsr(X, R, shifts=None, wX=None, wR=None, keep=None, thresh=1e-12): Weights applied by TSR. """ - from psutil import Process ndims = X.ndim X = unsqueeze(X) R = unsqueeze(R) @@ -184,7 +183,7 @@ def tsr(X, R, shifts=None, wX=None, wR=None, keep=None, thresh=1e-12): for t in np.arange(n_trials_X): r = multishift(R[..., t], shifts, reshape=True) y[..., t] = X[:z.shape[0], :, t] - (r @ regression) - + y, mean2 = demean(y, wX, return_mean=True, inplace=True) idx = np.arange(offset1, initial_samples - offset2) mean_total = mean1 + mean2 From c998c8131b25875411690d04b0c11d7fe76ab6ae Mon Sep 17 00:00:00 2001 From: Eduard Ort Date: Thu, 2 Dec 2021 16:22:45 +0100 Subject: [PATCH 3/6] fix pep and flake --- meegkit/dss.py | 6 +++--- meegkit/utils/denoise.py | 2 +- meegkit/utils/matrix.py | 2 +- meegkit/utils/sig.py | 4 ++-- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/meegkit/dss.py b/meegkit/dss.py index 87e915c5..2a1a8a3b 100644 --- a/meegkit/dss.py +++ b/meegkit/dss.py @@ -201,7 +201,7 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, # Cancel line_frequency and harmonics + light lowpass X_filt = smooth(X, sfreq / fline) - + # X - X_filt results in the artifact plus some residual biological signal # Reduce dimensionality to avoid overfitting if nkeep is not None: @@ -211,7 +211,7 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, else: X_noise_pca = (X - X_filt).copy() nkeep = n_chans - + # Compute blockwise covariances of raw and biased data n_harm = np.floor((sfreq / 2) / fline).astype(int) c0 = np.zeros((nkeep, nkeep)) @@ -225,7 +225,7 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, # bias data c0 += tscov(X_block)[0] c1 += tscov(gaussfilt(X_block, sfreq, fline, fwhm=1, n_harm=n_harm))[0] - + # DSS to isolate line components from residual todss, _, pwr0, pwr1 = dss0(c0, c1) diff --git a/meegkit/utils/denoise.py b/meegkit/utils/denoise.py index 9181397b..24ca2092 100644 --- a/meegkit/utils/denoise.py +++ b/meegkit/utils/denoise.py @@ -17,7 +17,7 @@ def demean(X, weights=None, return_mean=False, inplace=False): weights : array, shape=(n_samples) return_mean : bool, optional inplace : bool, optional - Save the resulting array in X, or in a new array + Save the resulting array in X, or in a new array Returns ------- diff --git a/meegkit/utils/matrix.py b/meegkit/utils/matrix.py index c7aa85fe..153d29f9 100644 --- a/meegkit/utils/matrix.py +++ b/meegkit/utils/matrix.py @@ -504,7 +504,7 @@ def fold(X, epoch_size): n_chans = X.shape[0] // epoch_size if X.shape[0] / epoch_size >= 1: X = np.reshape(X, (epoch_size, X.shape[1], n_chans), order="F") - #X = np.transpose(np.reshape(X, (epoch_size, n_chans, X.shape[1]), + # X = np.transpose(np.reshape(X, (epoch_size, n_chans, X.shape[1]), # order="F").copy(), [0, 2, 1]) return X diff --git a/meegkit/utils/sig.py b/meegkit/utils/sig.py index 499ca272..ada58d5d 100644 --- a/meegkit/utils/sig.py +++ b/meegkit/utils/sig.py @@ -347,10 +347,10 @@ def gaussfilt(data, srate, f, fwhm, n_harm=1, shift=0, return_empvals=False, # filter if data.ndim == 2: filtdat = 2 * np.real(np.fft.ifft( - np.fft.fft(data, axis=0)*fx[:, None], axis=0)) + np.fft.fft(data, axis=0) * fx[:, None], axis=0)) elif data.ndim == 3: filtdat = 2 * np.real(np.fft.ifft( - np.fft.fft(data, axis=0)*fx[:, None, None], + np.fft.fft(data, axis=0) * fx[:, None, None], axis=0)) if return_empvals or show: From d5f29e28fd5818073a55ba5248cc49d0744d6d91 Mon Sep 17 00:00:00 2001 From: Eduard Ort Date: Thu, 2 Dec 2021 16:34:36 +0100 Subject: [PATCH 4/6] undo matrix change --- meegkit/utils/matrix.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/meegkit/utils/matrix.py b/meegkit/utils/matrix.py index 153d29f9..51ab2a7e 100644 --- a/meegkit/utils/matrix.py +++ b/meegkit/utils/matrix.py @@ -503,9 +503,9 @@ def fold(X, epoch_size): n_chans = X.shape[0] // epoch_size if X.shape[0] / epoch_size >= 1: - X = np.reshape(X, (epoch_size, X.shape[1], n_chans), order="F") - # X = np.transpose(np.reshape(X, (epoch_size, n_chans, X.shape[1]), - # order="F").copy(), [0, 2, 1]) + # X = np.reshape(X, (epoch_size, X.shape[1], n_chans), order="F") + X = np.transpose(np.reshape(X, (epoch_size, n_chans, X.shape[1]), + order="F").copy(), [0, 2, 1]) return X From 76b0e485a629682f56d27d9bdf706c652ba83784 Mon Sep 17 00:00:00 2001 From: Eduard Ort Date: Thu, 2 Dec 2021 16:53:15 +0100 Subject: [PATCH 5/6] undo matrix reshaping --- meegkit/utils/matrix.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/meegkit/utils/matrix.py b/meegkit/utils/matrix.py index 51ab2a7e..01c79d2b 100644 --- a/meegkit/utils/matrix.py +++ b/meegkit/utils/matrix.py @@ -503,10 +503,9 @@ def fold(X, epoch_size): n_chans = X.shape[0] // epoch_size if X.shape[0] / epoch_size >= 1: - # X = np.reshape(X, (epoch_size, X.shape[1], n_chans), order="F") X = np.transpose(np.reshape(X, (epoch_size, n_chans, X.shape[1]), order="F").copy(), [0, 2, 1]) - return X + return X, Y def unfold(X): From 0a46afcc3a3440c1373090bcab70698899b00d27 Mon Sep 17 00:00:00 2001 From: nbara <10333715+nbara@users.noreply.github.com> Date: Thu, 2 Dec 2021 17:43:10 +0100 Subject: [PATCH 6/6] compromise --- meegkit/detrend.py | 10 ++++++---- meegkit/dss.py | 13 ++++++++----- meegkit/utils/denoise.py | 28 ++++++++++++++-------------- meegkit/utils/matrix.py | 11 ++++++----- tests/test_detrend.py | 2 +- tests/test_dss.py | 8 +++----- 6 files changed, 38 insertions(+), 34 deletions(-) diff --git a/meegkit/detrend.py b/meegkit/detrend.py index 0b882e25..34de9f74 100644 --- a/meegkit/detrend.py +++ b/meegkit/detrend.py @@ -200,7 +200,7 @@ def regress(x, r, w=None, threshold=1e-7, return_mean=False): xx = demean(x[:, i], wc) * wc # remove channel-specific-weighted mean from regressor - r = demean(r, wc) + r = demean(r, wc, inplace=True) rr = r * wc V, _ = pca(rr.T @ rr, thresh=threshold) rr = rr.dot(V) @@ -282,12 +282,14 @@ def _plot_detrend(x, y, w): f = plt.figure() gs = GridSpec(4, 1, figure=f) ax1 = f.add_subplot(gs[:3, 0]) - plt.plot(x, label='original', color='C0') - plt.plot(y, label='detrended', color='C1') + lines = ax1.plot(x, label='original', color='C0') + plt.setp(lines[1:], label="_") + lines = ax1.plot(y, label='detrended', color='C1') + plt.setp(lines[1:], label="_") ax1.set_xlim(0, n_times) ax1.set_xticklabels('') ax1.set_title('Robust detrending') - ax1.legend() + ax1.legend(fontsize='smaller') ax2 = f.add_subplot(gs[3, 0]) ax2.pcolormesh(w.T, cmap='Greys') diff --git a/meegkit/dss.py b/meegkit/dss.py index 2a1a8a3b..f4841ba3 100644 --- a/meegkit/dss.py +++ b/meegkit/dss.py @@ -192,7 +192,7 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, if X.shape[0] < nfft: print('Reducing nfft to {}'.format(X.shape[0])) nfft = X.shape[0] - n_samples, n_chans, n_trials = theshapeof(X) + n_samples, n_chans, _ = theshapeof(X) if blocksize is None: blocksize = n_samples @@ -203,13 +203,15 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, X_filt = smooth(X, sfreq / fline) # X - X_filt results in the artifact plus some residual biological signal + X_noise = X - X_filt + # Reduce dimensionality to avoid overfitting if nkeep is not None: - cov_X_res = tscov(X - X_filt)[0] + cov_X_res = tscov(X_noise)[0] V, _ = pca(cov_X_res, nkeep) - X_noise_pca = (X - X_filt) @ V + X_noise_pca = X_noise @ V else: - X_noise_pca = (X - X_filt).copy() + X_noise_pca = X_noise.copy() nkeep = n_chans # Compute blockwise covariances of raw and biased data @@ -240,7 +242,8 @@ def dss_line(X, fline, sfreq, nremove=1, nfft=1024, nkeep=None, blocksize=None, # Remove line components from X_noise idx_remove = np.arange(nremove) X_artifact = matmul3d(X_noise_pca, todss[:, idx_remove]) - X_res = tsr(X - X_filt, X_artifact)[0] # project them out + X_res = tsr(X_noise, X_artifact)[0] # project them out + # reconstruct clean signal y = X_filt + X_res diff --git a/meegkit/utils/denoise.py b/meegkit/utils/denoise.py index 24ca2092..92ed6321 100644 --- a/meegkit/utils/denoise.py +++ b/meegkit/utils/denoise.py @@ -15,13 +15,14 @@ def demean(X, weights=None, return_mean=False, inplace=False): X : array, shape=(n_samples, n_channels[, n_trials]) Data. weights : array, shape=(n_samples) - return_mean : bool, optional - inplace : bool, optional - Save the resulting array in X, or in a new array + return_mean : bool + If True, also return signal mean (default=False). + inplace : bool + If True, save the resulting array in X (default=False). Returns ------- - demeaned_X : array, shape=(n_samples, n_channels[, n_trials]) + X : array, shape=(n_samples, n_channels[, n_trials]) Centered data. mn : array Mean value. @@ -29,6 +30,9 @@ def demean(X, weights=None, return_mean=False, inplace=False): """ weights = _check_weights(weights, X) ndims = X.ndim + + if not inplace: + X = X.copy() n_samples, n_chans, n_trials = theshapeof(X) X = unfold(X) @@ -49,21 +53,17 @@ def demean(X, weights=None, return_mean=False, inplace=False): else: mn = np.mean(X, axis=0, keepdims=True) - # return new array or remove the mean in X itself - if inplace: - np.subtract(X, mn, out=X) - # create alias, so that rest of script can run with the same variable - demeaned_X = X - else: - demeaned_X = X - mn + # Remove mean (no copy) + X -= mn + # np.subtract(X, mn, out=X) if n_trials > 1 or ndims == 3: - demeaned_X = fold(demeaned_X, n_samples) + X = fold(X, n_samples) if return_mean: - return demeaned_X, mn # the_mean.shape=(1, the_mean.shape[0]) + return X, mn # the_mean.shape=(1, the_mean.shape[0]) else: - return demeaned_X + return X def mean_over_trials(X, weights=None): diff --git a/meegkit/utils/matrix.py b/meegkit/utils/matrix.py index 01c79d2b..e186b095 100644 --- a/meegkit/utils/matrix.py +++ b/meegkit/utils/matrix.py @@ -495,17 +495,18 @@ def unsqueeze(X): def fold(X, epoch_size): - """Fold 2D X into 3D.""" + """Fold 2D (n_times, n_channels) X into 3D (n_times, n_chans, n_trials).""" if X.ndim == 1: X = X[:, np.newaxis] if X.ndim > 2: raise AttributeError('X must be 2D at most') - n_chans = X.shape[0] // epoch_size + nt = X.shape[0] // epoch_size + nc = X.shape[1] if X.shape[0] / epoch_size >= 1: - X = np.transpose(np.reshape(X, (epoch_size, n_chans, X.shape[1]), - order="F").copy(), [0, 2, 1]) - return X, Y + return X.reshape((epoch_size, nt, nc), order="F").transpose([0, 2, 1]) + else: + return X def unfold(X): diff --git a/tests/test_detrend.py b/tests/test_detrend.py index 3419ec97..5d115044 100644 --- a/tests/test_detrend.py +++ b/tests/test_detrend.py @@ -11,7 +11,7 @@ def test_regress(): # Simple regression example, no weights # fit random walk y = np.cumsum(np.random.randn(1000, 1), axis=0) - x = np.arange(1000)[:, None] + x = np.arange(1000.)[:, None] x = np.hstack([x, x ** 2, x ** 3]) [b, z] = regress(y, x) diff --git a/tests/test_dss.py b/tests/test_dss.py index b15d794e..be390355 100644 --- a/tests/test_dss.py +++ b/tests/test_dss.py @@ -39,10 +39,10 @@ def test_dss0(n_bad_chans): atol=1e-6) # use abs as DSS component might be flipped -def test_dss1(show=False): +def test_dss1(show=True): """Test DSS1 (evoked).""" n_samples = 300 - data, source = create_line_data(n_samples=n_samples) + data, source = create_line_data(n_samples=n_samples, fline=.01) todss, _, pwr0, pwr1 = dss.dss1(data, weights=None, ) z = fold(np.dot(unfold(data), todss), epoch_size=n_samples) @@ -179,12 +179,10 @@ def profile_dss_line(nkeep): sr = 200 fline = 20 nsamples = 1000000 - nchans = 16 + nchans = 99 s = create_line_data(n_samples=3 * nsamples, n_chans=nchans, n_trials=1, fline=fline / sr, SNR=2)[0][..., 0] - # 2D case, n_outputs == 1 - pr = cProfile.Profile() pr.enable() out, _ = dss.dss_line(s, fline, sr, nkeep=nkeep)