From 4eda7759f8b1ab1f48b24fcc0843603b495b0378 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 26 Aug 2018 17:04:47 -0400 Subject: [PATCH 01/21] Add godec. --- tedana/decomposition/__init__.py | 5 +- tedana/decomposition/godecomp.py | 176 +++++++++++++++++++++++++++++++ tedana/workflows/tedana.py | 25 +++-- 3 files changed, 198 insertions(+), 8 deletions(-) create mode 100644 tedana/decomposition/godecomp.py diff --git a/tedana/decomposition/__init__.py b/tedana/decomposition/__init__.py index 72625ca0f..3dd102614 100644 --- a/tedana/decomposition/__init__.py +++ b/tedana/decomposition/__init__.py @@ -4,8 +4,11 @@ from .eigendecomp import ( tedpca, tedica, ) +from .godecomp import ( + tedgodec +) __all__ = [ - 'tedpca', 'tedica', + 'tedpca', 'tedica', 'tedgodec', ] diff --git a/tedana/decomposition/godecomp.py b/tedana/decomposition/godecomp.py new file mode 100644 index 000000000..e6f54f78e --- /dev/null +++ b/tedana/decomposition/godecomp.py @@ -0,0 +1,176 @@ +""" +Go Decomposition +""" +import logging + +import numpy as np +from numpy.linalg import qr, lstsq + +from tedana import utils +from tedana.decomposition._utils import dwtmat, idwtmat + +LGR = logging.getLogger(__name__) + + +def wthresh(a, thresh): + """ + Soft wavelet threshold + """ + res = np.abs(a) - thresh + return np.sign(a) * ((res > 0) * res) + + +def godec(data, thresh=.03, rank=2, power=1, tol=1e-3, max_iter=100, + random_seed=0, verbose=True): + """ + Perform Go Decomposition + + Default threshold of .03 is assumed to be for input in the range 0-1... + original matlab had 8 out of 255, which is about .03 scaled to 0-1 range + + Parameters + ---------- + data : (M x T) array_like + + Returns + ------- + L : array_like + Low-rank components. Similar to global signals. Should be discarded, + according to Power et al. (2018). + S : array_like + Sparse components. Should be retained, according to Power et al. + (2018). + G : array_like + Residuals (i.e., data minus sparse and low-rank components) + """ + LGR.info('Starting Go Decomposition') + _, n_vols = data.shape + L = data + S = np.zeros(L.shape) + itr = 0 + random_state = np.random.RandomState(random_seed) + while True: + Y2 = random_state.randn(n_vols, rank) + for i in range(power+1): + Y1 = np.dot(L, Y2) + Y2 = np.dot(L.T, Y1) + Q, R = qr(Y2) + L_new = np.dot(np.dot(L, Q), Q.T) + T = L - L_new + S + L = L_new + S = wthresh(T, thresh) + T -= S + err = np.linalg.norm(T.ravel(), 2) + if err < tol: + if verbose: + LGR.info('Successful convergence after {0} ' + 'iterations'.format(itr+1)) + break + elif itr >= max_iter: + if verbose: + LGR.warning('Model failed to converge after {0} ' + 'iterations'.format(itr+1)) + break + L += T + itr += 1 + + # Is this even useful in soft GoDec? May be a display issue... + G = data - L - S + return L, S, G + + +def _tedgodec(data, wavelet=False, rank=2, increment=2, power=2, tol=1e-3, + thresh=10, max_iter=500, norm_mode='vn', random_seed=0, + verbose=True): + """ + Perform TE-dependent Go Decomposition + + Parameters + ---------- + data : (M x T) array_like + wavelet : :obj:`bool`, optional + + """ + if norm_mode == 'dm': + # Demean + data_mean = data.mean(-1) + data_norm = data - data_mean[:, np.newaxis] + elif norm_mode == 'vn': + # What is this? + data_mean = data.mean(-1)[:, np.newaxis] + data_std = data.std(-1)[:, np.newaxis] + data_norm = (data - data_mean) / data_std + else: + data_norm = data + + # GoDec + if wavelet: + data_wt, cal = dwtmat(data_norm) + L, S, G = godec(data_wt, thresh=data_wt.std()*thresh, rank=rank, + power=power, tol=tol, max_iter=max_iter, + random_seed=random_seed, verbose=verbose) + L = idwtmat(L, cal) + S = idwtmat(S, cal) + G = idwtmat(G, cal) + else: + L, S, G = godec(data_norm, thresh=thresh, rank=rank, + power=power, tol=tol, max_iter=max_iter, + random_seed=random_seed, verbose=verbose) + + if norm_mode == 'dm': + # Remean + L += data_mean + elif norm_mode == 'vn': + L = (L * data_std) + data_mean + S *= data_std + G *= data_std + + return L, S, G + + +def tedgodec(optcom_ts, mmix, mask, acc, ref_img, ranks=[2], wavelet=False, + thresh=10, norm_mode='vn', power=2): + """ + optcom_ts : (S x T) array_like + Optimally combined time series data + mmix : (C x T) array_like + Mixing matrix for converting input data to component space, where `C` + is components and `T` is the same as in `optcom_ts` + mask : (S,) array_like + Boolean mask array + acc : :obj:`list` + Indices of accepted (BOLD) components in `mmix` + ref_img : :obj:`str` or img_like + Reference image to dictate how outputs are saved to disk + ranks : list of int + Ranks of low-rank components to run + """ + # Construct denoised data from optcom, mmix, and acc + optcom_masked = optcom_ts[mask, :] + optcom_mu = optcom_masked.mean(axis=-1)[:, np.newaxis] + optcom_std = optcom_masked.std(axis=-1)[:, np.newaxis] + data_norm = (optcom_masked - optcom_mu) / optcom_std + cbetas = lstsq(mmix, data_norm.T, rcond=None)[0].T + resid = data_norm - np.dot(cbetas, mmix.T) + bold_ts = np.dot(cbetas[:, acc], mmix[:, acc].T) + medn_ts = optcom_mu + ((bold_ts + resid) * optcom_std) + + for rank in ranks: + L, S, G = _tedgodec(medn_ts, rank=rank, power=power, thresh=thresh, + max_iter=500, norm_mode=norm_mode) + + if norm_mode is None: + name_norm_mode = '' + else: + name_norm_mode = 'n{0}'.format(norm_mode) + + if wavelet: + name_norm_mode = 'w{0}'.format(name_norm_mode) + + suffix = '{0}r{1}p{2}t{3}'.format(name_norm_mode, rank, power, thresh) + utils.filewrite(utils.unmask(L, mask), + 'lowrank_{0}.nii'.format(suffix), ref_img) + utils.filewrite(utils.unmask(S, mask), + 'sparse_{0}.nii'.format(suffix), ref_img) + utils.filewrite(utils.unmask(G, mask), + 'noise_{0}.nii'.format(suffix), ref_img) diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index dcd63b0c6..dab4167da 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -170,8 +170,9 @@ def _get_parser(): def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, - strict=False, gscontrol=True, kdaw=10., rdaw=1., conv=2.5e-5, - ste=-1, combmode='t2s', dne=False, + strict=False, gscontrol=True, ws_denoise=None, + kdaw=10., rdaw=1., + conv=2.5e-5, ste=-1, combmode='t2s', dne=False, initcost='tanh', finalcost='tanh', stabilize=False, filecsdata=False, wvpca=False, label=None, fixed_seed=42, debug=False, quiet=False): @@ -201,6 +202,8 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, Ignore low-variance ambiguous components. Default is False. gscontrol : :obj:`bool`, optional Control global signal using spatial approach. Default is True. + ws_denoise : {None, 'gsr', 'godec'}, optional + Which method to apply for widespread signal denoising. Default is None. kdaw : :obj:`float`, optional Dimensionality augmentation weight (Kappa). Default is 10. -1 for low-dimensional ICA. @@ -384,14 +387,14 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, utils.filewrite(s0G, op.join(out_dir, 's0vG.nii'), ref_img) # optimally combine data - OCcatd = model.make_optcom(catd, tes, mask, t2s=t2sG, combmode=combmode) + optcom_ts = model.make_optcom(catd, tes, mask, t2s=t2sG, combmode=combmode) # regress out global signal unless explicitly not desired if gscontrol: - catd, OCcatd = model.gscontrol_raw(catd, OCcatd, n_echos, ref_img) + catd, optcom_ts = model.gscontrol_raw(catd, optcom_ts, n_echos, ref_img) if mixm is None: - n_components, dd = decomposition.tedpca(catd, OCcatd, combmode, mask, + n_components, dd = decomposition.tedpca(catd, optcom_ts, combmode, mask, t2s, t2sG, stabilize, ref_img, tes=tes, kdaw=kdaw, rdaw=rdaw, ste=ste, wvpca=wvpca) @@ -430,9 +433,17 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, LGR.warning('No BOLD components detected! Please check data and ' 'results!') - utils.writeresults(OCcatd, mask, comptable, mmix, fixed_seed, n_vols, + utils.writeresults(optcom_ts, mask, comptable, mmix, fixed_seed, n_vols, acc, rej, midk, empty, ref_img) - utils.gscontrol_mmix(OCcatd, mmix, mask, acc, ref_img) + + # Widespread noise control + if ws_denoise == 'gsr': + utils.gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img) + elif ws_denoise == 'godec': + decomposition.tedgodec(optcom_ts, mmix, mask, acc, ref_img, + ranks=[2], wavelet=wvpca, + thresh=10, norm_mode='vn', power=2) + if dne: utils.writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img) From 0f6eb2fab38c68f98e5649f681a143d9b8b2d744 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 26 Aug 2018 18:15:31 -0400 Subject: [PATCH 02/21] Fix integration test workflow call. --- tedana/tests/test_tedana.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tedana/tests/test_tedana.py b/tedana/tests/test_tedana.py index 4a21e1eeb..979df6a8b 100644 --- a/tedana/tests/test_tedana.py +++ b/tedana/tests/test_tedana.py @@ -17,7 +17,7 @@ def test_basic_tedana(): files. """ workflows.tedana_workflow([op.expanduser('~/data/zcat_ffd.nii.gz')], - [14.5, 38.5, 62.5]) + [14.5, 38.5, 62.5], ws_denoise='gsr') assert op.isfile('comp_table.txt') From 6c4abc33a017b196b957515e1831af54bc408f4a Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Mon, 27 Aug 2018 22:49:29 -0400 Subject: [PATCH 03/21] Fix potential bug in T1c and GODEC. MEDN time series was being constructed without the ignored components. Also the High-Kappa time series in T1c was being constructed without the mean. --- tedana/decomposition/godecomp.py | 38 ++++++++++++++--------- tedana/model/fit.py | 2 +- tedana/utils/io.py | 53 +++++++++++++++++++------------- tedana/workflows/tedana.py | 4 +-- 4 files changed, 58 insertions(+), 39 deletions(-) diff --git a/tedana/decomposition/godecomp.py b/tedana/decomposition/godecomp.py index e6f54f78e..f239991fe 100644 --- a/tedana/decomposition/godecomp.py +++ b/tedana/decomposition/godecomp.py @@ -2,6 +2,7 @@ Go Decomposition """ import logging +import os.path as op import numpy as np from numpy.linalg import qr, lstsq @@ -12,7 +13,7 @@ LGR = logging.getLogger(__name__) -def wthresh(a, thresh): +def _wthresh(a, thresh): """ Soft wavelet threshold """ @@ -58,18 +59,17 @@ def godec(data, thresh=.03, rank=2, power=1, tol=1e-3, max_iter=100, L_new = np.dot(np.dot(L, Q), Q.T) T = L - L_new + S L = L_new - S = wthresh(T, thresh) + S = _wthresh(T, thresh) T -= S err = np.linalg.norm(T.ravel(), 2) if err < tol: if verbose: - LGR.info('Successful convergence after {0} ' - 'iterations'.format(itr+1)) + LGR.info('Successful convergence after %i iterations', itr+1) break elif itr >= max_iter: if verbose: - LGR.warning('Model failed to converge after {0} ' - 'iterations'.format(itr+1)) + LGR.warning('Model failed to converge after %i iterations', + itr+1) break L += T itr += 1 @@ -79,7 +79,7 @@ def godec(data, thresh=.03, rank=2, power=1, tol=1e-3, max_iter=100, return L, S, G -def _tedgodec(data, wavelet=False, rank=2, increment=2, power=2, tol=1e-3, +def _tedgodec(data, wavelet=False, rank=2, power=2, tol=1e-3, thresh=10, max_iter=500, norm_mode='vn', random_seed=0, verbose=True): """ @@ -96,7 +96,7 @@ def _tedgodec(data, wavelet=False, rank=2, increment=2, power=2, tol=1e-3, data_mean = data.mean(-1) data_norm = data - data_mean[:, np.newaxis] elif norm_mode == 'vn': - # What is this? + # Variance normalize data_mean = data.mean(-1)[:, np.newaxis] data_std = data.std(-1)[:, np.newaxis] data_norm = (data - data_mean) / data_std @@ -128,8 +128,8 @@ def _tedgodec(data, wavelet=False, rank=2, increment=2, power=2, tol=1e-3, return L, S, G -def tedgodec(optcom_ts, mmix, mask, acc, ref_img, ranks=[2], wavelet=False, - thresh=10, norm_mode='vn', power=2): +def tedgodec(optcom_ts, mmix, mask, acc, ign, ref_img, ranks=[2], + wavelet=False, thresh=10, norm_mode='vn', power=2, out_dir='.'): """ optcom_ts : (S x T) array_like Optimally combined time series data @@ -140,18 +140,23 @@ def tedgodec(optcom_ts, mmix, mask, acc, ref_img, ranks=[2], wavelet=False, Boolean mask array acc : :obj:`list` Indices of accepted (BOLD) components in `mmix` + ign : :obj:`list` + Indices of all ignored components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk ranks : list of int Ranks of low-rank components to run + norm_mode : {'vn', 'dm', None} """ - # Construct denoised data from optcom, mmix, and acc + # Construct denoised data from optcom, mmix, acc, and all_ref optcom_masked = optcom_ts[mask, :] optcom_mu = optcom_masked.mean(axis=-1)[:, np.newaxis] optcom_std = optcom_masked.std(axis=-1)[:, np.newaxis] data_norm = (optcom_masked - optcom_mu) / optcom_std cbetas = lstsq(mmix, data_norm.T, rcond=None)[0].T - resid = data_norm - np.dot(cbetas, mmix.T) + all_comps = np.arange(mmix.shape[0]) + not_ign = sorted(np.setdiff1d(all_comps, ign)) + resid = data_norm - np.dot(cbetas[:, not_ign], mmix[:, not_ign].T) bold_ts = np.dot(cbetas[:, acc], mmix[:, acc].T) medn_ts = optcom_mu + ((bold_ts + resid) * optcom_std) @@ -169,8 +174,11 @@ def tedgodec(optcom_ts, mmix, mask, acc, ref_img, ranks=[2], wavelet=False, suffix = '{0}r{1}p{2}t{3}'.format(name_norm_mode, rank, power, thresh) utils.filewrite(utils.unmask(L, mask), - 'lowrank_{0}.nii'.format(suffix), ref_img) + op.join(out_dir, 'lowrank_{0}.nii'.format(suffix)), + ref_img) utils.filewrite(utils.unmask(S, mask), - 'sparse_{0}.nii'.format(suffix), ref_img) + op.join(out_dir, 'sparse_{0}.nii'.format(suffix)), + ref_img) utils.filewrite(utils.unmask(G, mask), - 'noise_{0}.nii'.format(suffix), ref_img) + op.join(out_dir, 'noise_{0}.nii'.format(suffix)), + ref_img) diff --git a/tedana/model/fit.py b/tedana/model/fit.py index c7036a9d9..462bce507 100644 --- a/tedana/model/fit.py +++ b/tedana/model/fit.py @@ -376,7 +376,7 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, dtrank=4): Removes global signal from individual echo `catd` and `optcom` time series This function uses the spatial global signal estimation approach to - to removal global signal out of individual echo time series datasets. The + to remove global signal out of individual echo time series datasets. The spatial global signal is estimated from the optimally combined data after detrending with a Legendre polynomial basis of `order = 0` and `degree = dtrank`. diff --git a/tedana/utils/io.py b/tedana/utils/io.py index 333594a80..f1913dce9 100644 --- a/tedana/utils/io.py +++ b/tedana/utils/io.py @@ -13,7 +13,7 @@ LGR = logging.getLogger(__name__) -def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img): +def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): """ Perform global signal regression. @@ -28,6 +28,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img): Boolean mask array acc : :obj:`list` Indices of accepted (BOLD) components in `mmix` + ign : :obj:`list` + Indices of all ignored components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk @@ -39,7 +41,7 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img): Filename Content ====================== ================================================= sphis_hik.nii T1-like effect. - hik_ts_OC_T1c.nii T1 corrected time series. + hik_ts_OC_T1c.nii T1 corrected BOLD (high-Kappa) time series. dn_ts_OC_T1c.nii Denoised version of T1 corrected time series betas_hik_OC_T1c.nii T1-GS corrected components meica_mix_T1c.1D T1-GS corrected mixing matrix @@ -54,7 +56,9 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img): """ data_norm = (optcom_masked - optcom_mu) / optcom_std cbetas = lstsq(mmix, data_norm.T, rcond=None)[0].T - resid = data_norm - np.dot(cbetas, mmix.T) + all_comps = np.arange(mmix.shape[0]) + not_ign = sorted(np.setdiff1d(all_comps, ign)) + resid = data_norm - np.dot(cbetas[:, not_ign], mmix[:, not_ign].T) """ Build BOLD time series without amplitudes, and save T1-like effect @@ -75,14 +79,14 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img): """ bold_noT1gs = bold_ts - np.dot(lstsq(glob_sig.T, bold_ts.T, rcond=None)[0].T, glob_sig) - utils.filewrite(utils.unmask(bold_noT1gs * optcom_std, mask), - 'hik_ts_OC_T1c.nii', ref_img) + hik_ts = optcom_mu + (bold_noT1gs * optcom_std) + utils.filewrite(utils.unmask(hik_ts, mask), 'hik_ts_OC_T1c.nii', ref_img) """ Make denoised version of T1-corrected time series """ medn_ts = optcom_mu + ((bold_noT1gs + resid) * optcom_std) - utils.filewrite(utils.unmask(medn_ts, mask), 'dn_ts_OC_T1c', ref_img) + utils.filewrite(utils.unmask(medn_ts, mask), 'dn_ts_OC_T1c.nii', ref_img) """ Orthogonalize mixing matrix w.r.t. T1-GS @@ -99,8 +103,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img): Write T1-GS corrected components and mixing matrix """ cbetas_norm = lstsq(mmixnogs_norm.T, data_norm.T, rcond=None)[0].T - utils.filewrite(utils.unmask(cbetas_norm[:, 2:], mask), 'betas_hik_OC_T1c', - ref_img) + utils.filewrite(utils.unmask(cbetas_norm[:, 2:], mask), + 'betas_hik_OC_T1c.nii', ref_img) np.savetxt('meica_mix_T1c.1D', mmixnogs) @@ -177,10 +181,13 @@ def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix=''): ====================== ================================================= Filename Content ====================== ================================================= - hik_ts_[suffix].nii High-Kappa time series. + hik_ts_[suffix].nii High-Kappa time series. Combination of accepted + components. midk_ts_[suffix].nii Mid-Kappa time series. - low_ts_[suffix].nii Low-Kappa time series. - dn_ts_[suffix].nii Denoised time series. + low_ts_[suffix].nii Low-Kappa time series. Combination of rejected + components. + dn_ts_[suffix].nii Denoised time series. Optimally combined data + minus mid-Kappa and rejected components. ====================== ================================================= """ @@ -264,8 +271,8 @@ def writefeats(data, mmix, mask, ref_img, suffix=''): return fname -def writect(comptable, n_vols, fixed_seed, acc, rej, midk, empty, ctname='comp_table.txt', - varexpl='-1'): +def writect(comptable, n_vols, fixed_seed, acc, rej, midk, ign, + ctname='comp_table.txt', varexpl='-1'): """ Saves component table to disk @@ -285,7 +292,7 @@ def writect(comptable, n_vols, fixed_seed, acc, rej, midk, empty, ctname='comp_t Indices of rejected (non-BOLD) components in `mmix` midk : :obj:`list` Indices of mid-K (questionable) components in `mmix` - empty : :obj:`list` + ign : :obj:`list` Indices of ignored components in `mmix` ctname : :obj:`str`, optional Filename to save comptable to disk. Default 'comp_table.txt' @@ -303,8 +310,10 @@ def writect(comptable, n_vols, fixed_seed, acc, rej, midk, empty, ctname='comp_t indices. rejected.txt A comma-separated list of the rejected component indices. - midk_rejected.txt A comma-separated list of middle-kappa components + midk_rejected.txt A comma-separated list of middle-Kappa components which were ultimately rejected. + ignored.txt A comma-separated list of the ignored component + indices. [ctname] (comp_table.txt) Component table file. ========================= ================================================= """ @@ -321,6 +330,9 @@ def writect(comptable, n_vols, fixed_seed, acc, rej, midk, empty, ctname='comp_t with open('midk_rejected.txt', 'w') as fo: fo.write(','.join([str(int(cc)) for cc in midk])) + with open('ignored.txt', 'w') as fo: + fo.write(','.join([str(int(cc)) for cc in ign])) + _computed_vars = dict(file=op.abspath(op.curdir), vex=varexpl, seed=fixed_seed, @@ -331,7 +343,7 @@ def writect(comptable, n_vols, fixed_seed, acc, rej, midk, empty, ctname='comp_t acc=','.join([str(int(cc)) for cc in acc]), rej=','.join([str(int(cc)) for cc in rej]), mid=','.join([str(int(cc)) for cc in midk]), - ign=','.join([str(int(cc)) for cc in empty])) + ign=','.join([str(int(cc)) for cc in ign])) heading = textwrap.dedent("""\ # ME-ICA Component statistics table for: {file} # # Dataset variance explained by ICA (VEx): {vex:.2f} @@ -360,7 +372,7 @@ def writect(comptable, n_vols, fixed_seed, acc, rej, midk, empty, ctname='comp_t def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, - acc, rej, midk, empty, ref_img): + acc, rej, midk, ign, ref_img): """ Denoises `ts` and saves all resulting files to disk @@ -387,7 +399,7 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, Indices of rejected (non-BOLD) components in `mmix` midk : :obj:`list` Indices of mid-K (questionable) components in `mmix` - empty : :obj:`list` + ign : :obj:`list` Indices of ignored components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk @@ -433,8 +445,8 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, mmix[:, acc], mask, ref_img, suffix='OC2') LGR.info('Writing Z-normalized spatial component maps: {}'.format(op.abspath(fout))) - writect(comptable, n_vols, fixed_seed, acc, rej, midk, empty, ctname='comp_table.txt', - varexpl=varexpl) + writect(comptable, n_vols, fixed_seed, acc, rej, midk, ign, + ctname='comp_table.txt', varexpl=varexpl) LGR.info('Writing component table: {}'.format(op.abspath('comp_table.txt'))) @@ -481,7 +493,6 @@ def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img): :py:func:`tedana.utils.io.write_split_ts`. ====================== ================================================= """ - for i_echo in range(catd.shape[1]): LGR.info('Writing Kappa-filtered echo #{:01d} timeseries'.format(i_echo+1)) write_split_ts(catd[:, i_echo, :], mmix, mask, acc, rej, midk, ref_img, diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index dab4167da..4ae615d37 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -438,9 +438,9 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, # Widespread noise control if ws_denoise == 'gsr': - utils.gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img) + utils.gscontrol_mmix(optcom_ts, mmix, mask, acc, empty, ref_img) elif ws_denoise == 'godec': - decomposition.tedgodec(optcom_ts, mmix, mask, acc, ref_img, + decomposition.tedgodec(optcom_ts, mmix, mask, acc, empty, ref_img, ranks=[2], wavelet=wvpca, thresh=10, norm_mode='vn', power=2) From 60ce2f39c463eb9723d435d85762632395b4eea4 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Mon, 27 Aug 2018 23:00:34 -0400 Subject: [PATCH 04/21] Remove mean from High-Kappa time series. --- tedana/utils/io.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tedana/utils/io.py b/tedana/utils/io.py index f1913dce9..25de82928 100644 --- a/tedana/utils/io.py +++ b/tedana/utils/io.py @@ -79,7 +79,7 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): """ bold_noT1gs = bold_ts - np.dot(lstsq(glob_sig.T, bold_ts.T, rcond=None)[0].T, glob_sig) - hik_ts = optcom_mu + (bold_noT1gs * optcom_std) + hik_ts = bold_noT1gs * optcom_std utils.filewrite(utils.unmask(hik_ts, mask), 'hik_ts_OC_T1c.nii', ref_img) """ From b7119eaf39dacba3f5967f0d9682a59a9b264709 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 28 Aug 2018 07:07:33 -0400 Subject: [PATCH 05/21] Fix mmix shape bug. --- tedana/utils/io.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tedana/utils/io.py b/tedana/utils/io.py index 25de82928..bb056568c 100644 --- a/tedana/utils/io.py +++ b/tedana/utils/io.py @@ -21,7 +21,7 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): ---------- optcom_ts : (S x T) array_like Optimally combined time series data - mmix : (C x T) array_like + mmix : (T x C) array_like Mixing matrix for converting input data to component space, where `C` is components and `T` is the same as in `optcom_ts` mask : (S,) array_like @@ -56,7 +56,7 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): """ data_norm = (optcom_masked - optcom_mu) / optcom_std cbetas = lstsq(mmix, data_norm.T, rcond=None)[0].T - all_comps = np.arange(mmix.shape[0]) + all_comps = np.arange(mmix.shape[1]) not_ign = sorted(np.setdiff1d(all_comps, ign)) resid = data_norm - np.dot(cbetas[:, not_ign], mmix[:, not_ign].T) From dba5c24a19328dc8e9bbb4ba7dbcf3c549d61c40 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 28 Aug 2018 09:21:05 -0400 Subject: [PATCH 06/21] Improve gscontrol_mmix docstring. --- tedana/utils/io.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/tedana/utils/io.py b/tedana/utils/io.py index bb056568c..da0d32785 100644 --- a/tedana/utils/io.py +++ b/tedana/utils/io.py @@ -15,7 +15,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): """ - Perform global signal regression. + Perform minimum image regression to identify and remove global artifacts. + This is an alternative to methods like GoDec with fewer parameters to tune. Parameters ---------- @@ -40,11 +41,11 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): ====================== ================================================= Filename Content ====================== ================================================= - sphis_hik.nii T1-like effect. - hik_ts_OC_T1c.nii T1 corrected BOLD (high-Kappa) time series. - dn_ts_OC_T1c.nii Denoised version of T1 corrected time series - betas_hik_OC_T1c.nii T1-GS corrected components - meica_mix_T1c.1D T1-GS corrected mixing matrix + sphis_hik.nii T1-like effect + hik_ts_OC_T1c.nii T1-corrected BOLD (high-Kappa) time series + dn_ts_OC_T1c.nii Denoised version of T1-corrected time series + betas_hik_OC_T1c.nii T1 global signal-corrected components + meica_mix_T1c.1D T1 global signal-corrected mixing matrix ====================== ================================================= """ optcom_masked = optcom_ts[mask, :] From b1851a0952baa38d825cb7304ba91c644527d7d3 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 30 Aug 2018 21:31:28 -0400 Subject: [PATCH 07/21] Incorporate "out_dir" variable into workflow and functions. Still outputs to TED though. --- tedana/decomposition/eigendecomp.py | 17 +++-- tedana/model/fit.py | 27 ++++++-- tedana/selection/select_comps.py | 14 ++-- tedana/utils/io.py | 103 +++++++++++++++++++--------- tedana/workflows/tedana.py | 28 +++++--- 5 files changed, 130 insertions(+), 59 deletions(-) diff --git a/tedana/decomposition/eigendecomp.py b/tedana/decomposition/eigendecomp.py index 4b47ad3fa..bdcf060b2 100644 --- a/tedana/decomposition/eigendecomp.py +++ b/tedana/decomposition/eigendecomp.py @@ -19,7 +19,8 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, - ref_img, tes, kdaw, rdaw, ste=0, mlepca=True, wvpca=False): + ref_img, tes, kdaw, rdaw, ste=0, mlepca=True, wvpca=False, + out_dir='.'): """ Use principal components analysis (PCA) to identify and remove thermal noise from multi-echo data. @@ -57,6 +58,8 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, guessing PCA dimensionality instead of a traditional SVD. Default: True wvpca : :obj:`bool`, optional Whether to apply wavelet denoising to data. Default: False + out_dir : :obj:`str`, optional + Output directory in which to save output files Returns ------- @@ -137,7 +140,7 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, if wvpca: dz, cAl = dwtmat(dz) - if not op.exists('pcastate.pkl'): + if not op.exists(op.join(out_dir, 'pcastate.pkl')): # do PC dimension selection and get eigenvalue cutoff if mlepca: from sklearn.decomposition import PCA @@ -177,7 +180,7 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, ctb = np.vstack([ctb.T[:3], sp]).T # Save state - fname = op.abspath('pcastate.pkl') + fname = op.abspath(op.join(out_dir, 'pcastate.pkl')) LGR.info('Saving PCA results to: {}'.format(fname)) pcastate = {'u': u, 's': s, 'v': v, 'ctb': ctb, 'eigelb': eigelb, 'spmin': spmin, 'spcum': spcum} @@ -189,14 +192,16 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, else: # if loading existing state LGR.info('Loading PCA from: pcastate.pkl') - with open('pcastate.pkl', 'rb') as handle: + with open(op.join(out_dir, 'pcastate.pkl'), 'rb') as handle: pcastate = pickle.load(handle) u, s, v = pcastate['u'], pcastate['s'], pcastate['v'] ctb, eigelb = pcastate['ctb'], pcastate['eigelb'] spmin, spcum = pcastate['spmin'], pcastate['spcum'] - np.savetxt('comp_table_pca.txt', ctb[ctb[:, 1].argsort(), :][::-1]) - np.savetxt('mepca_mix.1D', v[ctb[:, 1].argsort()[::-1], :].T) + np.savetxt(op.join(out_dir, 'comp_table_pca.txt'), + ctb[ctb[:, 1].argsort(), :][::-1]) + np.savetxt(op.join(out_dir, 'mepca_mix.1D'), + v[ctb[:, 1].argsort()[::-1], :].T) kappas = ctb[ctb[:, 1].argsort(), 1] rhos = ctb[ctb[:, 2].argsort(), 2] diff --git a/tedana/model/fit.py b/tedana/model/fit.py index 462bce507..14b033ff8 100644 --- a/tedana/model/fit.py +++ b/tedana/model/fit.py @@ -2,6 +2,7 @@ Fit models. """ import logging +import os.path as op import nilearn.image as niimg from nilearn._utils import check_niimg @@ -371,7 +372,7 @@ def get_coeffs(data, X, mask=None, add_const=False): return betas -def gscontrol_raw(catd, optcom, n_echos, ref_img, dtrank=4): +def gscontrol_raw(catd, optcom, n_echos, ref_img, out_dir='.', dtrank=4): """ Removes global signal from individual echo `catd` and `optcom` time series @@ -391,6 +392,8 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, dtrank=4): Number of echos in data. Should be the same as `E` dimension of `catd` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk + out_dir : :obj:`str` + Output directory in which to save output files dtrank : :obj:`int`, optional Specifies degree of Legendre polynomial basis function for estimating spatial global signal. Default: 4 @@ -401,6 +404,19 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, dtrank=4): Input `catd` with global signal removed from time series dm_optcom : (S x T) array_like Input `optcom` with global signal removed from time series + + Notes + ----- + This function writes out several files. Files are listed below: + + ====================== ================================================= + Filename Content + ====================== ================================================= + T1gs.nii Spatial global signal + tsoc_orig.nii Optimally combined data + tsoc_nogs.nii Optimally combined data with global signal + regressed out + ====================== ================================================= """ LGR.info('Applying amplitude-based T1 equilibration correction') if catd.shape[0] != optcom.shape[0]: @@ -429,13 +445,14 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, dtrank=4): detr = dat - np.dot(sol.T, Lmix.T)[0] sphis = (detr).min(axis=1) sphis -= sphis.mean() - utils.filewrite(utils.unmask(sphis, Gmask), 'T1gs', ref_img) + utils.filewrite(utils.unmask(sphis, Gmask), + op.join(out_dir, 'T1gs.nii'), ref_img) # find time course ofc the spatial global signal # make basis with the Legendre basis glsig = np.linalg.lstsq(np.atleast_2d(sphis).T, dat, rcond=None)[0] glsig = stats.zscore(glsig, axis=None) - np.savetxt('glsig.1D', glsig) + np.savetxt(op.join(out_dir, 'glsig.1D'), glsig) glbase = np.hstack([Lmix, glsig.T]) # Project global signal out of optimally combined data @@ -443,9 +460,9 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, dtrank=4): tsoc_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T, np.atleast_2d(glbase.T[dtrank])) + Gmu[Gmask][:, np.newaxis] - utils.filewrite(optcom, 'tsoc_orig', ref_img) + utils.filewrite(optcom, op.join(out_dir, 'tsoc_orig.nii'), ref_img) dm_optcom = utils.unmask(tsoc_nogs, Gmask) - utils.filewrite(dm_optcom, 'tsoc_nogs', ref_img) + utils.filewrite(dm_optcom, op.join(out_dir, 'tsoc_nogs.nii'), ref_img) # Project glbase out of each echo dm_catd = catd.copy() # don't overwrite catd diff --git a/tedana/selection/select_comps.py b/tedana/selection/select_comps.py index 26d7d690e..36384462b 100644 --- a/tedana/selection/select_comps.py +++ b/tedana/selection/select_comps.py @@ -2,6 +2,7 @@ Functions to identify TE-dependent and TE-independent components. """ import os +import os.path as op import json import logging import pickle @@ -21,7 +22,8 @@ def selcomps(seldict, mmix, mask, ref_img, manacc, n_echos, t2s, s0, olevel=2, - oversion=99, filecsdata=True, savecsdiag=True, strict_mode=False): + oversion=99, filecsdata=True, savecsdiag=True, strict_mode=False, + out_dir='.'): """ Labels ICA components to keep or remove from denoised data @@ -60,6 +62,8 @@ def selcomps(seldict, mmix, mask, ref_img, manacc, n_echos, t2s, s0, olevel=2, Default: True strict_mode: :obj:`bool`, optional Default: False + out_dir : :obj:`str` + Output directory in which to save output files Returns ------- @@ -776,7 +780,9 @@ def selcomps(seldict, mmix, mask, ref_img, manacc, n_echos, t2s, s0, olevel=2, group0_res = np.intersect1d(KRguess, group0) phys_var_zs.append((vvex - vvex[group0_res].mean()) / vvex[group0_res].std()) veinBout = utils.unmask(veinmaskB, mask) - utils.filewrite(veinBout.astype(float), 'veins_l%i' % t2sl_i, ref_img) + utils.filewrite(veinBout.astype(float), + op.join(out_dir, 'veins_l{}.nii'.format(t2sl_i)), + ref_img) # Mask to sample veins phys_var_z = np.array(phys_var_zs).max(0) @@ -855,10 +861,10 @@ def selcomps(seldict, mmix, mask, ref_img, manacc, n_echos, t2s, s0, olevel=2, list(field_art), list(phys_art), list(misc_art), list(acc_comps), list(ign)] - with open('csstepdata.json', 'w') as ofh: + with open(op.join(out_dir, 'csstepdata.json'), 'w') as ofh: json.dump(dict(zip(diagstep_keys, diagstep_vals)), ofh, indent=4, sort_keys=True, default=str) allfz = np.array([Tz, Vz, Ktz, KRr, cnz, Rz, mmix_kurt, fdist_z]) - np.savetxt('csdata.txt', allfz) + np.savetxt(op.join(out_dir, 'csdata.txt'), allfz) return list(sorted(acc_comps)), list(sorted(rej)), list(sorted(midk)), list(sorted(ign)) diff --git a/tedana/utils/io.py b/tedana/utils/io.py index da0d32785..020d786b3 100644 --- a/tedana/utils/io.py +++ b/tedana/utils/io.py @@ -2,8 +2,8 @@ Functions to handle file input/output """ import logging -import os.path as op import textwrap +import os.path as op import numpy as np from numpy.linalg import lstsq @@ -13,7 +13,7 @@ LGR = logging.getLogger(__name__) -def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): +def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img, out_dir='.'): """ Perform minimum image regression to identify and remove global artifacts. This is an alternative to methods like GoDec with fewer parameters to tune. @@ -33,6 +33,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): Indices of all ignored components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk + out_dir : :obj:`str` + Output directory in which to save output files Notes ----- @@ -67,7 +69,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): bold_ts = np.dot(cbetas[:, acc], mmix[:, acc].T) t1_map = bold_ts.min(axis=-1) t1_map -= t1_map.mean() - utils.filewrite(utils.unmask(t1_map, mask), 'sphis_hik', ref_img) + utils.filewrite(utils.unmask(t1_map, mask), + op.join(out_dir, 'sphis_hik.nii'), ref_img) t1_map = t1_map[:, np.newaxis] """ @@ -81,13 +84,15 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): bold_noT1gs = bold_ts - np.dot(lstsq(glob_sig.T, bold_ts.T, rcond=None)[0].T, glob_sig) hik_ts = bold_noT1gs * optcom_std - utils.filewrite(utils.unmask(hik_ts, mask), 'hik_ts_OC_T1c.nii', ref_img) + utils.filewrite(utils.unmask(hik_ts, mask), + op.join(out_dir, 'hik_ts_OC_T1c.nii'), ref_img) """ Make denoised version of T1-corrected time series """ medn_ts = optcom_mu + ((bold_noT1gs + resid) * optcom_std) - utils.filewrite(utils.unmask(medn_ts, mask), 'dn_ts_OC_T1c.nii', ref_img) + utils.filewrite(utils.unmask(medn_ts, mask), + op.join(out_dir, 'dn_ts_OC_T1c.nii'), ref_img) """ Orthogonalize mixing matrix w.r.t. T1-GS @@ -105,8 +110,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img): """ cbetas_norm = lstsq(mmixnogs_norm.T, data_norm.T, rcond=None)[0].T utils.filewrite(utils.unmask(cbetas_norm[:, 2:], mask), - 'betas_hik_OC_T1c.nii', ref_img) - np.savetxt('meica_mix_T1c.1D', mmixnogs) + op.join(out_dir, 'betas_hik_OC_T1c.nii'), ref_img) + np.savetxt(op.join(out_dir, 'meica_mix_T1c.1D'), mmixnogs) def split_ts(data, mmix, mask, acc): @@ -146,7 +151,8 @@ def split_ts(data, mmix, mask, acc): return hikts, resid -def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix=''): +def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix='', + out_dir='.'): """ Splits `data` into denoised / noise / ignored time series and saves to disk @@ -169,6 +175,8 @@ def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix=''): Reference image to dictate how outputs are saved to disk suffix : :obj:`str`, optional Appended to name of saved files (before extension). Default: '' + out_dir : :obj:`str` + Output directory in which to save output files Returns ------- @@ -211,27 +219,34 @@ def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix=''): if len(acc) != 0: fout = utils.filewrite(utils.unmask(hikts, mask), - 'hik_ts_{0}'.format(suffix), ref_img) + op.join(out_dir, + 'hik_ts_{0}.nii'.format(suffix)), + ref_img) LGR.info('Writing high-Kappa time series: {}'.format(op.abspath(fout))) if len(midk) != 0: fout = utils.filewrite(utils.unmask(midkts, mask), - 'midk_ts_{0}'.format(suffix), ref_img) + op.join(out_dir, + 'midk_ts_{0}.nii'.format(suffix)), + ref_img) LGR.info('Writing mid-Kappa time series: {}'.format(op.abspath(fout))) if len(rej) != 0: fout = utils.filewrite(utils.unmask(lowkts, mask), - 'lowk_ts_{0}'.format(suffix), ref_img) + op.join(out_dir, + 'lowk_ts_{0}.nii'.format(suffix)), + ref_img) LGR.info('Writing low-Kappa time series: {}'.format(op.abspath(fout))) fout = utils.filewrite(utils.unmask(dnts, mask), - 'dn_ts_{0}'.format(suffix), ref_img) + op.join(out_dir, + 'dn_ts_{0}.nii'.format(suffix)), ref_img) LGR.info('Writing denoised time series: {}'.format(op.abspath(fout))) return varexpl -def writefeats(data, mmix, mask, ref_img, suffix=''): +def writefeats(data, mmix, mask, ref_img, suffix='', out_dir='.'): """ Converts `data` to component space with `mmix` and saves to disk @@ -248,6 +263,8 @@ def writefeats(data, mmix, mask, ref_img, suffix=''): Reference image to dictate how outputs are saved to disk suffix : :obj:`str`, optional Appended to name of saved files (before extension). Default: '' + out_dir : :obj:`str` + Output directory in which to save output files Returns ------- @@ -267,13 +284,15 @@ def writefeats(data, mmix, mask, ref_img, suffix=''): # write feature versions of components feats = utils.unmask(model.computefeats2(data, mmix, mask), mask) - fname = utils.filewrite(feats, 'feats_{0}'.format(suffix), ref_img) + fname = utils.filewrite(feats, op.join(out_dir, + 'feats_{0}.nii'.format(suffix)), + ref_img) return fname def writect(comptable, n_vols, fixed_seed, acc, rej, midk, ign, - ctname='comp_table.txt', varexpl='-1'): + ctname='comp_table.txt', varexpl='-1', out_dir='.'): """ Saves component table to disk @@ -299,6 +318,8 @@ def writect(comptable, n_vols, fixed_seed, acc, rej, midk, ign, Filename to save comptable to disk. Default 'comp_table.txt' varexpl : :obj:`str` Variance explained by original data + out_dir : :obj:`str` + Output directory in which to save output files Notes ----- @@ -322,16 +343,19 @@ def writect(comptable, n_vols, fixed_seed, acc, rej, midk, ign, n_components = comptable.shape[0] sortab = comptable[comptable[:, 1].argsort()[::-1], :] - with open('accepted.txt', 'w') as fo: + if not op.dirname(ctname): + ctname = op.join(out_dir, ctname) + + with open(op.join(out_dir, 'accepted.txt'), 'w') as fo: fo.write(','.join([str(int(cc)) for cc in acc])) - with open('rejected.txt', 'w') as fo: + with open(op.join(out_dir, 'rejected.txt'), 'w') as fo: fo.write(','.join([str(int(cc)) for cc in rej])) - with open('midk_rejected.txt', 'w') as fo: + with open(op.join(out_dir, 'midk_rejected.txt'), 'w') as fo: fo.write(','.join([str(int(cc)) for cc in midk])) - with open('ignored.txt', 'w') as fo: + with open(op.join(out_dir, 'ignored.txt'), 'w') as fo: fo.write(','.join([str(int(cc)) for cc in ign])) _computed_vars = dict(file=op.abspath(op.curdir), @@ -373,7 +397,7 @@ def writect(comptable, n_vols, fixed_seed, acc, rej, midk, ign, def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, - acc, rej, midk, ign, ref_img): + acc, rej, midk, ign, ref_img, out_dir='.'): """ Denoises `ts` and saves all resulting files to disk @@ -404,6 +428,8 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, Indices of ignored components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk + out_dir : :obj:`str` + Output directory in which to save output files Notes ----- @@ -430,28 +456,37 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, ====================== ================================================= """ - fout = utils.filewrite(ts, 'ts_OC', ref_img) - LGR.info('Writing optimally-combined time series: {}'.format(op.abspath(fout))) + fout = utils.filewrite(ts, op.join(out_dir, 'ts_OC.nii'), ref_img) + LGR.info('Writing optimally-combined time series: ' + '{}'.format(op.abspath(fout))) - varexpl = write_split_ts(ts, mmix, mask, acc, rej, midk, ref_img, suffix='OC') + varexpl = write_split_ts(ts, mmix, mask, acc, rej, midk, ref_img, + suffix='OC', out_dir=out_dir) ts_B = model.get_coeffs(ts, mmix, mask) - fout = utils.filewrite(ts_B, 'betas_OC', ref_img) - LGR.info('Writing full ICA coefficient feature set: {}'.format(op.abspath(fout))) + fout = utils.filewrite(ts_B, op.join(out_dir, 'betas_OC.nii'), ref_img) + LGR.info('Writing full ICA coefficient feature set: ' + '{}'.format(op.abspath(fout))) if len(acc) != 0: - fout = utils.filewrite(ts_B[:, acc], 'betas_hik_OC', ref_img) - LGR.info('Writing denoised ICA coefficient feature set: {}'.format(op.abspath(fout))) + fout = utils.filewrite(ts_B[:, acc], + op.join(out_dir, 'betas_hik_OC.nii'), ref_img) + LGR.info('Writing denoised ICA coefficient feature set: ' + '{}'.format(op.abspath(fout))) fout = writefeats(split_ts(ts, mmix, mask, acc)[0], - mmix[:, acc], mask, ref_img, suffix='OC2') - LGR.info('Writing Z-normalized spatial component maps: {}'.format(op.abspath(fout))) + mmix[:, acc], mask, ref_img, suffix='OC2', + out_dir=out_dir) + LGR.info('Writing Z-normalized spatial component maps: ' + '{}'.format(op.abspath(fout))) + ct_file = op.join(out_dir, 'comp_table.txt') writect(comptable, n_vols, fixed_seed, acc, rej, midk, ign, - ctname='comp_table.txt', varexpl=varexpl) - LGR.info('Writing component table: {}'.format(op.abspath('comp_table.txt'))) + ctname=ct_file, varexpl=varexpl, out_dir=out_dir) + LGR.info('Writing component table: {}'.format(op.abspath(ct_file))) -def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img): +def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img, + out_dir='.'): """ Saves individually denoised echos to disk @@ -472,6 +507,8 @@ def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img): Indices of mid-K (questionable) components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk + out_dir : :obj:`str` + Output directory in which to save output files Notes ----- @@ -497,7 +534,7 @@ def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img): for i_echo in range(catd.shape[1]): LGR.info('Writing Kappa-filtered echo #{:01d} timeseries'.format(i_echo+1)) write_split_ts(catd[:, i_echo, :], mmix, mask, acc, rej, midk, ref_img, - suffix='e%i' % (i_echo+1)) + suffix='e{}'.format(i_echo+1), out_dir=out_dir) def ctabsel(ctabfile): diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index 4ae615d37..96c730987 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -359,8 +359,6 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, elif ctab is not None: raise IOError('Argument "ctab" must be an existing file.') - os.chdir(out_dir) - if mask is None: LGR.info('Computing adaptive mask') else: @@ -391,13 +389,15 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, # regress out global signal unless explicitly not desired if gscontrol: - catd, optcom_ts = model.gscontrol_raw(catd, optcom_ts, n_echos, ref_img) + catd, optcom_ts = model.gscontrol_raw(catd, optcom_ts, n_echos, + ref_img, out_dir=out_dir) if mixm is None: n_components, dd = decomposition.tedpca(catd, optcom_ts, combmode, mask, t2s, t2sG, stabilize, ref_img, tes=tes, kdaw=kdaw, rdaw=rdaw, - ste=ste, wvpca=wvpca) + ste=ste, wvpca=wvpca, + out_dir=out_dir) mmix_orig, fixed_seed = decomposition.tedica(n_components, dd, conv, fixed_seed, cost=initcost, final_cost=finalcost, verbose=debug) @@ -410,9 +410,12 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, reindex=True) np.savetxt(op.join(out_dir, 'meica_mix.1D'), mmix) - acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask, ref_img, manacc, - n_echos, t2s, s0, strict_mode=strict, - filecsdata=filecsdata) + acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask, + ref_img, manacc, + n_echos, t2s, s0, + strict_mode=strict, + filecsdata=filecsdata, + out_dir=out_dir) else: LGR.info('Using supplied mixing matrix from ICA') mmix_orig = np.loadtxt(op.join(out_dir, 'meica_mix.1D')) @@ -434,18 +437,21 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, 'results!') utils.writeresults(optcom_ts, mask, comptable, mmix, fixed_seed, n_vols, - acc, rej, midk, empty, ref_img) + acc, rej, midk, empty, ref_img, out_dir=out_dir) # Widespread noise control if ws_denoise == 'gsr': - utils.gscontrol_mmix(optcom_ts, mmix, mask, acc, empty, ref_img) + utils.gscontrol_mmix(optcom_ts, mmix, mask, acc, empty, ref_img, + out_dir=out_dir) elif ws_denoise == 'godec': decomposition.tedgodec(optcom_ts, mmix, mask, acc, empty, ref_img, ranks=[2], wavelet=wvpca, - thresh=10, norm_mode='vn', power=2) + thresh=10, norm_mode='vn', power=2, + out_dir=out_dir) if dne: - utils.writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img) + utils.writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img, + out_dir=out_dir) def _main(argv=None): From e72af8f0542bbe256adc4bd0a04fdf13690282e4 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 31 Aug 2018 13:45:13 -0400 Subject: [PATCH 08/21] Fix test (ish) --- tedana/tests/test_tedana.py | 2 +- tedana/utils/io.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tedana/tests/test_tedana.py b/tedana/tests/test_tedana.py index 979df6a8b..958945535 100644 --- a/tedana/tests/test_tedana.py +++ b/tedana/tests/test_tedana.py @@ -18,7 +18,7 @@ def test_basic_tedana(): """ workflows.tedana_workflow([op.expanduser('~/data/zcat_ffd.nii.gz')], [14.5, 38.5, 62.5], ws_denoise='gsr') - assert op.isfile('comp_table.txt') + assert op.isfile('~/code/TED.zcat_ffd/comp_table.txt') def compare_nifti(fn, test_dir, res_dir): diff --git a/tedana/utils/io.py b/tedana/utils/io.py index 020d786b3..bb0d85c34 100644 --- a/tedana/utils/io.py +++ b/tedana/utils/io.py @@ -457,7 +457,7 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, """ fout = utils.filewrite(ts, op.join(out_dir, 'ts_OC.nii'), ref_img) - LGR.info('Writing optimally-combined time series: ' + LGR.info('Writing optimally combined time series: ' '{}'.format(op.abspath(fout))) varexpl = write_split_ts(ts, mmix, mask, acc, rej, midk, ref_img, From a63a00ec2f93ca4927147a9711ab9e0acc4b6e16 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 31 Aug 2018 17:30:16 -0400 Subject: [PATCH 09/21] Now *really* fix tests. Well, one test. --- tedana/tests/test_tedana.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tedana/tests/test_tedana.py b/tedana/tests/test_tedana.py index 958945535..2b45f7225 100644 --- a/tedana/tests/test_tedana.py +++ b/tedana/tests/test_tedana.py @@ -18,7 +18,7 @@ def test_basic_tedana(): """ workflows.tedana_workflow([op.expanduser('~/data/zcat_ffd.nii.gz')], [14.5, 38.5, 62.5], ws_denoise='gsr') - assert op.isfile('~/code/TED.zcat_ffd/comp_table.txt') + assert op.isfile(op.expanduser('~/code/TED.zcat_ffd/comp_table.txt')) def compare_nifti(fn, test_dir, res_dir): From b27dff3e1b054e21dab2f80530d01d2c6e9d2d9b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 31 Aug 2018 17:42:54 -0400 Subject: [PATCH 10/21] Fix t2smap and associated tests. Update CLI. --- tedana/tests/test_t2smap.py | 18 +++++++++--------- tedana/tests/test_tedana.py | 3 ++- tedana/workflows/t2smap.py | 32 +++++++++++--------------------- tedana/workflows/tedana.py | 31 ++++++++++--------------------- 4 files changed, 32 insertions(+), 52 deletions(-) diff --git a/tedana/tests/test_t2smap.py b/tedana/tests/test_t2smap.py index 4b24189fe..aefb2239e 100644 --- a/tedana/tests/test_t2smap.py +++ b/tedana/tests/test_t2smap.py @@ -18,12 +18,12 @@ def test_basic_t2smap1(self): files. """ data_dir = get_test_data_path() + out_dir = 'TED.t2smap' data = [op.join(data_dir, 'echo1.nii.gz'), op.join(data_dir, 'echo2.nii.gz'), op.join(data_dir, 'echo3.nii.gz')] workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='t2s', - fitmode='all', label='t2smap') - out_dir = 'TED.echo1.t2smap' + fitmode='all', out_dir=out_dir) # Check outputs assert op.isfile(op.join(out_dir, 'ts_OC.nii')) @@ -44,12 +44,12 @@ def test_basic_t2smap2(self): files when fitmode is set to ts. """ data_dir = get_test_data_path() + out_dir = 'TED.t2smap' data = [op.join(data_dir, 'echo1.nii.gz'), op.join(data_dir, 'echo2.nii.gz'), op.join(data_dir, 'echo3.nii.gz')] workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='t2s', - fitmode='ts', label='t2smap') - out_dir = 'TED.echo1.t2smap' + fitmode='ts', out_dir=out_dir) # Check outputs assert op.isfile(op.join(out_dir, 'ts_OC.nii')) @@ -70,12 +70,12 @@ def test_basic_t2smap3(self): files when combmode is set to 'ste'. """ data_dir = get_test_data_path() + out_dir = 'TED.t2smap' data = [op.join(data_dir, 'echo1.nii.gz'), op.join(data_dir, 'echo2.nii.gz'), op.join(data_dir, 'echo3.nii.gz')] workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='ste', - fitmode='all', label='t2smap') - out_dir = 'TED.echo1.t2smap' + fitmode='all', out_dir=out_dir) # Check outputs assert op.isfile(op.join(out_dir, 'ts_OC.nii')) @@ -98,12 +98,12 @@ def test_basic_t2smap4(self): Not sure why this fails. """ data_dir = get_test_data_path() + out_dir = 'TED.t2smap' data = [op.join(data_dir, 'echo1.nii.gz'), op.join(data_dir, 'echo2.nii.gz'), op.join(data_dir, 'echo3.nii.gz')] workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='ste', - fitmode='ts', label='t2smap') - out_dir = 'TED.echo1.t2smap' + fitmode='ts', out_dir=out_dir) # Check outputs assert op.isfile(op.join(out_dir, 'ts_OC.nii')) @@ -120,4 +120,4 @@ def test_basic_t2smap4(self): def teardown_method(self): # Clean up folders - rmtree('TED.echo1.t2smap') + rmtree('TED.t2smap') diff --git a/tedana/tests/test_tedana.py b/tedana/tests/test_tedana.py index 2b45f7225..0457d22b3 100644 --- a/tedana/tests/test_tedana.py +++ b/tedana/tests/test_tedana.py @@ -17,7 +17,8 @@ def test_basic_tedana(): files. """ workflows.tedana_workflow([op.expanduser('~/data/zcat_ffd.nii.gz')], - [14.5, 38.5, 62.5], ws_denoise='gsr') + [14.5, 38.5, 62.5], ws_denoise='gsr', + out_dir='~/code/TED.zcat_ffd') assert op.isfile(op.expanduser('~/code/TED.zcat_ffd/comp_table.txt')) diff --git a/tedana/workflows/t2smap.py b/tedana/workflows/t2smap.py index aa8c6277e..609fd94ff 100644 --- a/tedana/workflows/t2smap.py +++ b/tedana/workflows/t2smap.py @@ -67,15 +67,16 @@ def _get_parser(): help=('Combination scheme for TEs: ' 't2s (Posse 1999, default), ste (Poser)'), default='t2s') - parser.add_argument('--label', - dest='label', + parser.add_argument('--out', + dest='out_dir', type=str, - help='Label for output directory.', - default=None) + help='Output directory.', + default='.') return parser -def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s', label=None): +def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s', + out_dir='.'): """ Estimate T2 and S0, and optimally combine data across TEs. @@ -96,15 +97,13 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s', label=N Default is 'all'. combmode : {'t2s', 'ste'}, optional Combination scheme for TEs: 't2s' (Posse 1999, default), 'ste' (Poser). - label : :obj:`str` or :obj:`None`, optional - Label for output directory. Default is None. + out_dir : :obj:`str`, optional + Output directory in which to save output files. Default is current + directory ('.'). Notes ----- - This workflow writes out several files, which are written out to a folder - named TED.[ref_label].[label] if ``label`` is provided and TED.[ref_label] - if not. ``ref_label`` is determined based on the name of the first ``data`` - file. + This workflow writes out several files to ``out_dir``. Files are listed below: @@ -137,16 +136,7 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s', label=N n_samp, n_echos, n_vols = catd.shape LGR.debug('Resulting data shape: {}'.format(catd.shape)) - try: - ref_label = os.path.basename(ref_img).split('.')[0] - except TypeError: - ref_label = os.path.basename(str(data[0])).split('.')[0] - - if label is not None: - out_dir = 'TED.{0}.{1}'.format(ref_label, label) - else: - out_dir = 'TED.{0}'.format(ref_label) - out_dir = op.abspath(out_dir) + out_dir = op.abspath(op.expanduser(out_dir)) if not op.isdir(out_dir): LGR.info('Creating output directory: {}'.format(out_dir)) os.mkdir(out_dir) diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index 96c730987..924307df6 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -144,11 +144,11 @@ def _get_parser(): help='Perform PCA on wavelet-transformed data', action='store_true', default=False) - parser.add_argument('--label', - dest='label', + parser.add_argument('--out', + dest='out_dir', type=str, - help='Label for output directory.', - default=None) + help='Output directory.', + default='.') parser.add_argument('--seed', dest='fixed_seed', type=int, @@ -175,7 +175,7 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, conv=2.5e-5, ste=-1, combmode='t2s', dne=False, initcost='tanh', finalcost='tanh', stabilize=False, filecsdata=False, wvpca=False, - label=None, fixed_seed=42, debug=False, quiet=False): + out_dir='.', fixed_seed=42, debug=False, quiet=False): """ Run the "canonical" TE-Dependent ANAlysis workflow. @@ -231,8 +231,9 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, wvpca : :obj:`bool`, optional Whether or not to perform PCA on wavelet-transformed data. Default is False. - label : :obj:`str` or :obj:`None`, optional - Label for output directory. Default is None. + out_dir : :obj:`str`, optional + Output directory in which to save output files. Default is current + directory ('.'). Other Parameters ---------------- @@ -260,10 +261,7 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, PROCEDURE 2a: Model fitting and component selection routines - This workflow writes out several files, which are written out to a folder - named TED.[ref_label].[label] if ``label`` is provided and TED.[ref_label] - if not. ``ref_label`` is determined based on the name of the first ``data`` - file. + This workflow writes out several files to ``out_dir``. Files are listed below: @@ -331,16 +329,7 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, kdaw, rdaw = float(kdaw), float(rdaw) - try: - ref_label = op.basename(ref_img).split('.')[0] - except TypeError: - ref_label = op.basename(str(data[0])).split('.')[0] - - if label is not None: - out_dir = 'TED.{0}.{1}'.format(ref_label, label) - else: - out_dir = 'TED.{0}'.format(ref_label) - out_dir = op.abspath(out_dir) + out_dir = op.abspath(op.expanduser(out_dir)) if not op.isdir(out_dir): LGR.info('Creating output directory: {}'.format(out_dir)) os.mkdir(out_dir) From 21b5d9843f376470dc825eda9f52bc80bf3ca135 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 31 Aug 2018 18:23:50 -0400 Subject: [PATCH 11/21] Refactor tedana workflow. Make more readable. --- docs/api.rst | 1 + tedana/workflows/tedana.py | 70 +++++++++++++++++++++----------------- 2 files changed, 39 insertions(+), 32 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index 572ac28f3..f1d62c1da 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -60,6 +60,7 @@ API tedana.decomposition.tedpca tedana.decomposition.tedica + tedana.decomposition.tedgodec :template: module.rst diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index 924307df6..efd8e7b19 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -9,7 +9,9 @@ import argparse import numpy as np from scipy import stats -from tedana import (decomposition, model, selection, utils) +import tedana.decomposition as decomp +import tedana.selection as sel +from tedana import (model, utils) from tedana.workflows.parser_utils import is_valid_file LGR = logging.getLogger(__name__) @@ -358,23 +360,26 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, LGR.debug('Retaining {}/{} samples'.format(mask.sum(), n_samp)) LGR.info('Computing T2* map') - t2s, s0, t2ss, s0s, t2sG, s0G = model.fit_decay(catd, tes, mask, masksum) + (t2s_limited, s0_limited, + t2ss, s0s, + t2s_full, s0_full) = model.fit_decay(catd, tes, mask, masksum) # set a hard cap for the T2* map # anything that is 10x higher than the 99.5 %ile will be reset to 99.5 %ile - cap_t2s = stats.scoreatpercentile(t2s.flatten(), 99.5, + cap_t2s = stats.scoreatpercentile(t2s_limited.flatten(), 99.5, interpolation_method='lower') LGR.debug('Setting cap on T2* map at {:.5f}'.format(cap_t2s * 10)) - t2s[t2s > cap_t2s * 10] = cap_t2s - utils.filewrite(t2s, op.join(out_dir, 't2sv.nii'), ref_img) - utils.filewrite(s0, op.join(out_dir, 's0v.nii'), ref_img) + t2s_limited[t2s_limited > cap_t2s * 10] = cap_t2s + utils.filewrite(t2s_limited, op.join(out_dir, 't2sv.nii'), ref_img) + utils.filewrite(s0_limited, op.join(out_dir, 's0v.nii'), ref_img) utils.filewrite(t2ss, op.join(out_dir, 't2ss.nii'), ref_img) utils.filewrite(s0s, op.join(out_dir, 's0vs.nii'), ref_img) - utils.filewrite(t2sG, op.join(out_dir, 't2svG.nii'), ref_img) - utils.filewrite(s0G, op.join(out_dir, 's0vG.nii'), ref_img) + utils.filewrite(t2s_full, op.join(out_dir, 't2svG.nii'), ref_img) + utils.filewrite(s0_full, op.join(out_dir, 's0vG.nii'), ref_img) # optimally combine data - optcom_ts = model.make_optcom(catd, tes, mask, t2s=t2sG, combmode=combmode) + optcom_ts = model.make_optcom(catd, tes, mask, t2s=t2s_full, + combmode=combmode) # regress out global signal unless explicitly not desired if gscontrol: @@ -382,40 +387,41 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, ref_img, out_dir=out_dir) if mixm is None: - n_components, dd = decomposition.tedpca(catd, optcom_ts, combmode, mask, - t2s, t2sG, stabilize, ref_img, - tes=tes, kdaw=kdaw, rdaw=rdaw, - ste=ste, wvpca=wvpca, - out_dir=out_dir) - mmix_orig, fixed_seed = decomposition.tedica(n_components, dd, conv, fixed_seed, - cost=initcost, final_cost=finalcost, - verbose=debug) + n_components, dd = decomp.tedpca(catd, optcom_ts, combmode, mask, + t2s_limited, t2s_full, stabilize, + ref_img, + tes=tes, kdaw=kdaw, rdaw=rdaw, + ste=ste, wvpca=wvpca, + out_dir=out_dir) + mmix_orig, fixed_seed = decomp.tedica(n_components, dd, conv, fixed_seed, + cost=initcost, final_cost=finalcost, + verbose=debug) np.savetxt(op.join(out_dir, '__meica_mix.1D'), mmix_orig) LGR.info('Making second component selection guess from ICA results') - seldict, comptable, betas, mmix = model.fitmodels_direct(catd, mmix_orig, - mask, t2s, t2sG, - tes, combmode, - ref_img, - reindex=True) + (seldict, comptable, + betas, mmix) = model.fitmodels_direct(catd, mmix_orig, + mask, t2s_limited, t2s_full, + tes, combmode, + ref_img, reindex=True) np.savetxt(op.join(out_dir, 'meica_mix.1D'), mmix) acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask, ref_img, manacc, - n_echos, t2s, s0, + n_echos, t2s_limited, s0_limited, strict_mode=strict, filecsdata=filecsdata, out_dir=out_dir) else: LGR.info('Using supplied mixing matrix from ICA') mmix_orig = np.loadtxt(op.join(out_dir, 'meica_mix.1D')) - seldict, comptable, betas, mmix = model.fitmodels_direct(catd, mmix_orig, - mask, t2s, t2sG, - tes, combmode, - ref_img) + (seldict, comptable, + betas, mmix) = model.fitmodels_direct(catd, mmix_orig, + mask, t2s_limited, t2s_full, + tes, combmode, ref_img) if ctab is None: acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask, ref_img, manacc, - n_echos, t2s, s0, + n_echos, t2s_limited, s0_limited, filecsdata=filecsdata, strict_mode=strict) else: @@ -433,10 +439,10 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, utils.gscontrol_mmix(optcom_ts, mmix, mask, acc, empty, ref_img, out_dir=out_dir) elif ws_denoise == 'godec': - decomposition.tedgodec(optcom_ts, mmix, mask, acc, empty, ref_img, - ranks=[2], wavelet=wvpca, - thresh=10, norm_mode='vn', power=2, - out_dir=out_dir) + decomp.tedgodec(optcom_ts, mmix, mask, acc, empty, ref_img, + ranks=[2], wavelet=wvpca, + thresh=10, norm_mode='vn', power=2, + out_dir=out_dir) if dne: utils.writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img, From 60da7f3b5df9657ef38dbcf96989c296d0517ac1 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 1 Sep 2018 10:56:12 -0400 Subject: [PATCH 12/21] Fix import issue. --- tedana/workflows/tedana.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index efd8e7b19..9c18e2a99 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -405,12 +405,12 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, ref_img, reindex=True) np.savetxt(op.join(out_dir, 'meica_mix.1D'), mmix) - acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask, - ref_img, manacc, - n_echos, t2s_limited, s0_limited, - strict_mode=strict, - filecsdata=filecsdata, - out_dir=out_dir) + acc, rej, midk, empty = sel.selcomps(seldict, mmix, mask, + ref_img, manacc, + n_echos, t2s_limited, s0_limited, + strict_mode=strict, + filecsdata=filecsdata, + out_dir=out_dir) else: LGR.info('Using supplied mixing matrix from ICA') mmix_orig = np.loadtxt(op.join(out_dir, 'meica_mix.1D')) @@ -419,11 +419,11 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, mask, t2s_limited, t2s_full, tes, combmode, ref_img) if ctab is None: - acc, rej, midk, empty = selection.selcomps(seldict, mmix, mask, - ref_img, manacc, - n_echos, t2s_limited, s0_limited, - filecsdata=filecsdata, - strict_mode=strict) + acc, rej, midk, empty = sel.selcomps(seldict, mmix, mask, + ref_img, manacc, + n_echos, t2s_limited, s0_limited, + filecsdata=filecsdata, + strict_mode=strict) else: acc, rej, midk, empty = utils.ctabsel(ctab) From c6cfc8835d601708e0461bf0e973c1f3a1d845ac Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 9 Sep 2018 13:19:19 -0400 Subject: [PATCH 13/21] Revert changes to gscontrol_mmix. --- tedana/utils/io.py | 32 +++++++++++++------------------- tedana/workflows/tedana.py | 2 +- 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/tedana/utils/io.py b/tedana/utils/io.py index bb0d85c34..3759031e1 100644 --- a/tedana/utils/io.py +++ b/tedana/utils/io.py @@ -13,7 +13,7 @@ LGR = logging.getLogger(__name__) -def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img, out_dir='.'): +def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img, out_dir='.'): """ Perform minimum image regression to identify and remove global artifacts. This is an alternative to methods like GoDec with fewer parameters to tune. @@ -22,15 +22,13 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img, out_dir='.'): ---------- optcom_ts : (S x T) array_like Optimally combined time series data - mmix : (T x C) array_like + mmix : (C x T) array_like Mixing matrix for converting input data to component space, where `C` is components and `T` is the same as in `optcom_ts` mask : (S,) array_like Boolean mask array acc : :obj:`list` Indices of accepted (BOLD) components in `mmix` - ign : :obj:`list` - Indices of all ignored components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk out_dir : :obj:`str` @@ -43,11 +41,11 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img, out_dir='.'): ====================== ================================================= Filename Content ====================== ================================================= - sphis_hik.nii T1-like effect - hik_ts_OC_T1c.nii T1-corrected BOLD (high-Kappa) time series - dn_ts_OC_T1c.nii Denoised version of T1-corrected time series - betas_hik_OC_T1c.nii T1 global signal-corrected components - meica_mix_T1c.1D T1 global signal-corrected mixing matrix + sphis_hik.nii T1-like effect. + hik_ts_OC_T1c.nii T1 corrected time series. + dn_ts_OC_T1c.nii Denoised version of T1 corrected time series + betas_hik_OC_T1c.nii T1-GS corrected components + meica_mix_T1c.1D T1-GS corrected mixing matrix ====================== ================================================= """ optcom_masked = optcom_ts[mask, :] @@ -59,9 +57,7 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img, out_dir='.'): """ data_norm = (optcom_masked - optcom_mu) / optcom_std cbetas = lstsq(mmix, data_norm.T, rcond=None)[0].T - all_comps = np.arange(mmix.shape[1]) - not_ign = sorted(np.setdiff1d(all_comps, ign)) - resid = data_norm - np.dot(cbetas[:, not_ign], mmix[:, not_ign].T) + resid = data_norm - np.dot(cbetas, mmix.T) """ Build BOLD time series without amplitudes, and save T1-like effect @@ -83,16 +79,14 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img, out_dir='.'): """ bold_noT1gs = bold_ts - np.dot(lstsq(glob_sig.T, bold_ts.T, rcond=None)[0].T, glob_sig) - hik_ts = bold_noT1gs * optcom_std - utils.filewrite(utils.unmask(hik_ts, mask), - op.join(out_dir, 'hik_ts_OC_T1c.nii'), ref_img) + utils.filewrite(utils.unmask(bold_noT1gs * optcom_std, mask), + 'hik_ts_OC_T1c.nii', ref_img) """ Make denoised version of T1-corrected time series """ medn_ts = optcom_mu + ((bold_noT1gs + resid) * optcom_std) - utils.filewrite(utils.unmask(medn_ts, mask), - op.join(out_dir, 'dn_ts_OC_T1c.nii'), ref_img) + utils.filewrite(utils.unmask(medn_ts, mask), 'dn_ts_OC_T1c', ref_img) """ Orthogonalize mixing matrix w.r.t. T1-GS @@ -109,8 +103,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ign, ref_img, out_dir='.'): Write T1-GS corrected components and mixing matrix """ cbetas_norm = lstsq(mmixnogs_norm.T, data_norm.T, rcond=None)[0].T - utils.filewrite(utils.unmask(cbetas_norm[:, 2:], mask), - op.join(out_dir, 'betas_hik_OC_T1c.nii'), ref_img) + utils.filewrite(utils.unmask(cbetas_norm[:, 2:], mask), 'betas_hik_OC_T1c', + ref_img) np.savetxt(op.join(out_dir, 'meica_mix_T1c.1D'), mmixnogs) diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index 9c18e2a99..09d459ca9 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -436,7 +436,7 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, # Widespread noise control if ws_denoise == 'gsr': - utils.gscontrol_mmix(optcom_ts, mmix, mask, acc, empty, ref_img, + utils.gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img, out_dir=out_dir) elif ws_denoise == 'godec': decomp.tedgodec(optcom_ts, mmix, mask, acc, empty, ref_img, From 9fe412323f8bacee8d199ec428a7582221175191 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 9 Sep 2018 13:21:59 -0400 Subject: [PATCH 14/21] And... fix output directory after reverting changes.. --- tedana/utils/io.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tedana/utils/io.py b/tedana/utils/io.py index 3759031e1..e726ea368 100644 --- a/tedana/utils/io.py +++ b/tedana/utils/io.py @@ -80,13 +80,14 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img, out_dir='.'): bold_noT1gs = bold_ts - np.dot(lstsq(glob_sig.T, bold_ts.T, rcond=None)[0].T, glob_sig) utils.filewrite(utils.unmask(bold_noT1gs * optcom_std, mask), - 'hik_ts_OC_T1c.nii', ref_img) + op.join(out_dir, 'hik_ts_OC_T1c.nii'), ref_img) """ Make denoised version of T1-corrected time series """ medn_ts = optcom_mu + ((bold_noT1gs + resid) * optcom_std) - utils.filewrite(utils.unmask(medn_ts, mask), 'dn_ts_OC_T1c', ref_img) + utils.filewrite(utils.unmask(medn_ts, mask), + op.join(out_dir, 'dn_ts_OC_T1c.nii'), ref_img) """ Orthogonalize mixing matrix w.r.t. T1-GS @@ -103,7 +104,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, acc, ref_img, out_dir='.'): Write T1-GS corrected components and mixing matrix """ cbetas_norm = lstsq(mmixnogs_norm.T, data_norm.T, rcond=None)[0].T - utils.filewrite(utils.unmask(cbetas_norm[:, 2:], mask), 'betas_hik_OC_T1c', + utils.filewrite(utils.unmask(cbetas_norm[:, 2:], mask), + op.join(out_dir, 'betas_hik_OC_T1c.nii'), ref_img) np.savetxt(op.join(out_dir, 'meica_mix_T1c.1D'), mmixnogs) From cf449659dd65178a2197418056a9d074b8878495 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 18 Sep 2018 11:46:38 -0400 Subject: [PATCH 15/21] Fix style issue. --- tedana/selection/select_comps.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tedana/selection/select_comps.py b/tedana/selection/select_comps.py index 423d3a94d..6bf4092fb 100644 --- a/tedana/selection/select_comps.py +++ b/tedana/selection/select_comps.py @@ -1,11 +1,10 @@ """ Functions to identify TE-dependent and TE-independent components. """ -import os -import os.path as op import json import logging import pickle +import os.path as op from nilearn._utils import check_niimg import numpy as np From b69ff09c032d88ed0343da44abd565127706a45a Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 28 Nov 2018 17:43:34 -0500 Subject: [PATCH 16/21] Undo changes to gitignore. --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index 1ea329bd2..85bbe91f0 100644 --- a/.gitignore +++ b/.gitignore @@ -101,3 +101,6 @@ ENV/ # mypy .mypy_cache/ + + # vscode + .vscode From d838a4af7e725c1d7486340fdfcfaea44f68a3f1 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 28 Nov 2018 17:43:57 -0500 Subject: [PATCH 17/21] Fix. --- .gitignore | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 85bbe91f0..645abefdb 100644 --- a/.gitignore +++ b/.gitignore @@ -102,5 +102,5 @@ ENV/ # mypy .mypy_cache/ - # vscode - .vscode +# vscode +.vscode From 92c2c8c521459ea71c4474551f0fe9ffb4c14f42 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 28 Nov 2018 18:15:59 -0500 Subject: [PATCH 18/21] Revert out_dir changes. --- tedana/decomposition/eigendecomp.py | 9 ++-- tedana/io.py | 71 ++++++++++------------------- tedana/model/fit.py | 12 ++--- tedana/tests/test_t2smap.py | 18 ++++---- tedana/workflows/t2smap.py | 18 ++++---- tedana/workflows/tedana.py | 27 +++++------ 6 files changed, 64 insertions(+), 91 deletions(-) diff --git a/tedana/decomposition/eigendecomp.py b/tedana/decomposition/eigendecomp.py index 41ff8fd4b..d2e415cbe 100644 --- a/tedana/decomposition/eigendecomp.py +++ b/tedana/decomposition/eigendecomp.py @@ -60,8 +60,7 @@ def run_svd(data): def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, - ref_img, tes, kdaw, rdaw, ste=0, wvpca=False, - out_dir='.'): + ref_img, tes, kdaw, rdaw, ste=0, wvpca=False): """ Use principal components analysis (PCA) to identify and remove thermal noise from multi-echo data. @@ -96,8 +95,6 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, a subset of echos. Default: 0 wvpca : :obj:`bool`, optional Whether to apply wavelet denoising to data. Default: False - out_dir : :obj:`str`, optional - Output directory in which to save output files Returns ------- @@ -214,7 +211,7 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, ct_df['normalized variance explained'] = varex_norm # Save state - fname = op.abspath(op.join(out_dir, 'pcastate.pkl')) + fname = op.abspath('pcastate.pkl') LGR.info('Saving PCA results to: {}'.format(fname)) pcastate = {'voxel_comp_weights': voxel_comp_weights, 'varex': varex, @@ -231,7 +228,7 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, else: # if loading existing state LGR.info('Loading PCA from: pcastate.pkl') - with open(op.join(out_dir, 'pcastate.pkl'), 'rb') as handle: + with open('pcastate.pkl', 'rb') as handle: pcastate = pickle.load(handle) voxel_comp_weights, varex = pcastate['voxel_comp_weights'], pcastate['varex'] comp_ts = pcastate['comp_ts'] diff --git a/tedana/io.py b/tedana/io.py index 920807d85..fbdc7e88b 100644 --- a/tedana/io.py +++ b/tedana/io.py @@ -16,7 +16,7 @@ LGR = logging.getLogger(__name__) -def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img, out_dir='.'): +def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img): """ Perform minimum image regression to identify and remove global artifacts. This is an alternative to methods like GoDec with fewer parameters to tune. @@ -34,8 +34,6 @@ def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img, out_dir='.'): Indices of accepted (BOLD) components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk - out_dir : :obj:`str` - Output directory in which to save output files Notes ----- @@ -69,8 +67,7 @@ def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img, out_dir='.'): bold_ts = np.dot(cbetas[:, acc], mmix[:, acc].T) t1_map = bold_ts.min(axis=-1) t1_map -= t1_map.mean() - filewrite(utils.unmask(t1_map, mask), - op.join(out_dir, 'sphis_hik.nii'), ref_img) + filewrite(utils.unmask(t1_map, mask), 'sphis_hik', ref_img) t1_map = t1_map[:, np.newaxis] """ @@ -84,14 +81,13 @@ def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img, out_dir='.'): bold_noT1gs = bold_ts - np.dot(lstsq(glob_sig.T, bold_ts.T, rcond=None)[0].T, glob_sig) filewrite(utils.unmask(bold_noT1gs * optcom_std, mask), - op.join(out_dir, 'hik_ts_OC_T1c.nii'), ref_img) + 'hik_ts_OC_T1c.nii', ref_img) """ Make denoised version of T1-corrected time series """ medn_ts = optcom_mu + ((bold_noT1gs + resid) * optcom_std) - filewrite(utils.unmask(medn_ts, mask), - op.join(out_dir, 'dn_ts_OC_T1c.nii'), ref_img) + filewrite(utils.unmask(medn_ts, mask), 'dn_ts_OC_T1c', ref_img) """ Orthogonalize mixing matrix w.r.t. T1-GS @@ -108,10 +104,9 @@ def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img, out_dir='.'): Write T1-GS corrected components and mixing matrix """ cbetas_norm = lstsq(mmixnogs_norm.T, data_norm.T, rcond=None)[0].T - filewrite(utils.unmask(cbetas_norm[:, 2:], mask), - op.join(out_dir, 'betas_hik_OC_T1c.nii'), - ref_img) - np.savetxt(op.join(out_dir, 'meica_mix_T1c.1D'), mmixnogs) + filewrite(utils.unmask(cbetas_norm[:, 2:], mask), 'betas_hik_OC_T1c', + ref_img) + np.savetxt('meica_mix_T1c.1D', mmixnogs) def split_ts(data, mmix, mask, acc): @@ -151,8 +146,7 @@ def split_ts(data, mmix, mask, acc): return hikts, resid -def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix='', - out_dir='.'): +def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix=''): """ Splits `data` into denoised / noise / ignored time series and saves to disk @@ -175,8 +169,6 @@ def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix='', Reference image to dictate how outputs are saved to disk suffix : :obj:`str`, optional Appended to name of saved files (before extension). Default: '' - out_dir : :obj:`str` - Output directory in which to save output files Returns ------- @@ -219,34 +211,27 @@ def write_split_ts(data, mmix, mask, acc, rej, midk, ref_img, suffix='', if len(acc) != 0: fout = filewrite(utils.unmask(hikts, mask), - op.join(out_dir, - 'hik_ts_{0}.nii'.format(suffix)), - ref_img) + 'hik_ts_{0}'.format(suffix), ref_img) LGR.info('Writing high-Kappa time series: {}'.format(op.abspath(fout))) if len(midk) != 0: fout = filewrite(utils.unmask(midkts, mask), - op.join(out_dir, - 'midk_ts_{0}.nii'.format(suffix)), - ref_img) + 'midk_ts_{0}'.format(suffix), ref_img) LGR.info('Writing mid-Kappa time series: {}'.format(op.abspath(fout))) if len(rej) != 0: fout = filewrite(utils.unmask(lowkts, mask), - op.join(out_dir, - 'lowk_ts_{0}.nii'.format(suffix)), - ref_img) + 'lowk_ts_{0}'.format(suffix), ref_img) LGR.info('Writing low-Kappa time series: {}'.format(op.abspath(fout))) fout = filewrite(utils.unmask(dnts, mask), - op.join(out_dir, - 'dn_ts_{0}.nii'.format(suffix)), ref_img) + 'dn_ts_{0}'.format(suffix), ref_img) LGR.info('Writing denoised time series: {}'.format(op.abspath(fout))) return varexpl -def writefeats(data, mmix, mask, ref_img, suffix='', out_dir='.'): +def writefeats(data, mmix, mask, ref_img, suffix=''): """ Converts `data` to component space with `mmix` and saves to disk @@ -263,8 +248,6 @@ def writefeats(data, mmix, mask, ref_img, suffix='', out_dir='.'): Reference image to dictate how outputs are saved to disk suffix : :obj:`str`, optional Appended to name of saved files (before extension). Default: '' - out_dir : :obj:`str` - Output directory in which to save output files Returns ------- @@ -284,15 +267,13 @@ def writefeats(data, mmix, mask, ref_img, suffix='', out_dir='.'): # write feature versions of components feats = utils.unmask(model.computefeats2(data, mmix, mask), mask) - fname = filewrite(feats, op.join(out_dir, - 'feats_{0}.nii'.format(suffix)), - ref_img) + fname = filewrite(feats, 'feats_{0}'.format(suffix), ref_img) return fname def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, - acc, rej, midk, ign, ref_img, out_dir='.'): + acc, rej, midk, ign, ref_img): """ Denoises `ts` and saves all resulting files to disk @@ -323,8 +304,6 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, Indices of ignored components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk - out_dir : :obj:`str` - Output directory in which to save output files Notes ----- @@ -351,31 +330,28 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, ====================== ================================================= """ - fout = filewrite(ts, op.join(out_dir, 'ts_OC.nii'), ref_img) + fout = filewrite(ts, 'ts_OC', ref_img) LGR.info('Writing optimally combined time series: ' '{}'.format(op.abspath(fout))) - varexpl = write_split_ts(ts, mmix, mask, acc, rej, midk, ref_img, - suffix='OC', out_dir=out_dir) + write_split_ts(ts, mmix, mask, acc, rej, midk, ref_img, suffix='OC') ts_B = model.get_coeffs(ts, mmix, mask) - fout = filewrite(ts_B, op.join(out_dir, 'betas_OC.nii'), ref_img) + fout = filewrite(ts_B, 'betas_OC', ref_img) LGR.info('Writing full ICA coefficient feature set: ' '{}'.format(op.abspath(fout))) if len(acc) != 0: - fout = filewrite(ts_B[:, acc], op.join(out_dir, 'betas_hik_OC.nii'), ref_img) + fout = filewrite(ts_B[:, acc], 'betas_hik_OC', ref_img) LGR.info('Writing denoised ICA coefficient feature set: ' '{}'.format(op.abspath(fout))) fout = writefeats(split_ts(ts, mmix, mask, acc)[0], - mmix[:, acc], mask, ref_img, suffix='OC2', - out_dir=out_dir) + mmix[:, acc], mask, ref_img, suffix='OC2') LGR.info('Writing Z-normalized spatial component maps: ' '{}'.format(op.abspath(fout))) -def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img, - out_dir='.'): +def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img): """ Saves individually denoised echos to disk @@ -396,8 +372,6 @@ def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img, Indices of mid-K (questionable) components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk - out_dir : :obj:`str` - Output directory in which to save output files Notes ----- @@ -420,10 +394,11 @@ def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img, :py:func:`tedana.utils.io.write_split_ts`. ====================== ================================================= """ + for i_echo in range(catd.shape[1]): LGR.info('Writing Kappa-filtered echo #{:01d} timeseries'.format(i_echo+1)) write_split_ts(catd[:, i_echo, :], mmix, mask, acc, rej, midk, ref_img, - suffix='e{}'.format(i_echo+1), out_dir=out_dir) + suffix='e{}'.format(i_echo+1)) def ctabsel(ctabfile): diff --git a/tedana/model/fit.py b/tedana/model/fit.py index 3bbb5798d..61b624db7 100644 --- a/tedana/model/fit.py +++ b/tedana/model/fit.py @@ -382,7 +382,7 @@ def get_coeffs(data, X, mask=None, add_const=False): return betas -def gscontrol_raw(catd, optcom, n_echos, ref_img, out_dir='.', dtrank=4): +def gscontrol_raw(catd, optcom, n_echos, ref_img, dtrank=4): """ Removes global signal from individual echo `catd` and `optcom` time series @@ -402,8 +402,6 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, out_dir='.', dtrank=4): Number of echos in data. Should be the same as `E` dimension of `catd` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk - out_dir : :obj:`str` - Output directory in which to save output files dtrank : :obj:`int`, optional Specifies degree of Legendre polynomial basis function for estimating spatial global signal. Default: 4 @@ -455,13 +453,13 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, out_dir='.', dtrank=4): detr = dat - np.dot(sol.T, Lmix.T)[0] sphis = (detr).min(axis=1) sphis -= sphis.mean() - io.filewrite(utils.unmask(sphis, Gmask), op.join(out_dir, 'T1gs.nii'), ref_img) + io.filewrite(utils.unmask(sphis, Gmask), 'T1gs', ref_img) # find time course ofc the spatial global signal # make basis with the Legendre basis glsig = np.linalg.lstsq(np.atleast_2d(sphis).T, dat, rcond=None)[0] glsig = stats.zscore(glsig, axis=None) - np.savetxt(op.join(out_dir, 'glsig.1D'), glsig) + np.savetxt('glsig.1D', glsig) glbase = np.hstack([Lmix, glsig.T]) # Project global signal out of optimally combined data @@ -469,9 +467,9 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, out_dir='.', dtrank=4): tsoc_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T, np.atleast_2d(glbase.T[dtrank])) + Gmu[Gmask][:, np.newaxis] - io.filewrite(optcom, op.join(out_dir, 'tsoc_orig.nii'), ref_img) + io.filewrite(optcom, 'tsoc_orig', ref_img) dm_optcom = utils.unmask(tsoc_nogs, Gmask) - io.filewrite(dm_optcom, op.join(out_dir, 'tsoc_nogs.nii'), ref_img) + io.filewrite(dm_optcom, 'tsoc_nogs', ref_img) # Project glbase out of each echo dm_catd = catd.copy() # don't overwrite catd diff --git a/tedana/tests/test_t2smap.py b/tedana/tests/test_t2smap.py index aefb2239e..4b24189fe 100644 --- a/tedana/tests/test_t2smap.py +++ b/tedana/tests/test_t2smap.py @@ -18,12 +18,12 @@ def test_basic_t2smap1(self): files. """ data_dir = get_test_data_path() - out_dir = 'TED.t2smap' data = [op.join(data_dir, 'echo1.nii.gz'), op.join(data_dir, 'echo2.nii.gz'), op.join(data_dir, 'echo3.nii.gz')] workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='t2s', - fitmode='all', out_dir=out_dir) + fitmode='all', label='t2smap') + out_dir = 'TED.echo1.t2smap' # Check outputs assert op.isfile(op.join(out_dir, 'ts_OC.nii')) @@ -44,12 +44,12 @@ def test_basic_t2smap2(self): files when fitmode is set to ts. """ data_dir = get_test_data_path() - out_dir = 'TED.t2smap' data = [op.join(data_dir, 'echo1.nii.gz'), op.join(data_dir, 'echo2.nii.gz'), op.join(data_dir, 'echo3.nii.gz')] workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='t2s', - fitmode='ts', out_dir=out_dir) + fitmode='ts', label='t2smap') + out_dir = 'TED.echo1.t2smap' # Check outputs assert op.isfile(op.join(out_dir, 'ts_OC.nii')) @@ -70,12 +70,12 @@ def test_basic_t2smap3(self): files when combmode is set to 'ste'. """ data_dir = get_test_data_path() - out_dir = 'TED.t2smap' data = [op.join(data_dir, 'echo1.nii.gz'), op.join(data_dir, 'echo2.nii.gz'), op.join(data_dir, 'echo3.nii.gz')] workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='ste', - fitmode='all', out_dir=out_dir) + fitmode='all', label='t2smap') + out_dir = 'TED.echo1.t2smap' # Check outputs assert op.isfile(op.join(out_dir, 'ts_OC.nii')) @@ -98,12 +98,12 @@ def test_basic_t2smap4(self): Not sure why this fails. """ data_dir = get_test_data_path() - out_dir = 'TED.t2smap' data = [op.join(data_dir, 'echo1.nii.gz'), op.join(data_dir, 'echo2.nii.gz'), op.join(data_dir, 'echo3.nii.gz')] workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='ste', - fitmode='ts', out_dir=out_dir) + fitmode='ts', label='t2smap') + out_dir = 'TED.echo1.t2smap' # Check outputs assert op.isfile(op.join(out_dir, 'ts_OC.nii')) @@ -120,4 +120,4 @@ def test_basic_t2smap4(self): def teardown_method(self): # Clean up folders - rmtree('TED.t2smap') + rmtree('TED.echo1.t2smap') diff --git a/tedana/workflows/t2smap.py b/tedana/workflows/t2smap.py index 6a113eb71..156719945 100644 --- a/tedana/workflows/t2smap.py +++ b/tedana/workflows/t2smap.py @@ -67,11 +67,11 @@ def _get_parser(): help=('Combination scheme for TEs: ' 't2s (Posse 1999, default), ste (Poser)'), default='t2s') - parser.add_argument('--out', - dest='out_dir', + parser.add_argument('--label', + dest='label', type=str, help='Label for output directory.', - default='.') + default=None) parser.add_argument('--debug', dest='debug', help=argparse.SUPPRESS, @@ -86,7 +86,7 @@ def _get_parser(): def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s', - out_dir='.', debug=False, quiet=False): + label=None, debug=False, quiet=False): """ Estimate T2 and S0, and optimally combine data across TEs. @@ -107,9 +107,8 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s', Default is 'all'. combmode : {'t2s', 'ste'}, optional Combination scheme for TEs: 't2s' (Posse 1999, default), 'ste' (Poser). - out_dir : :obj:`str`, optional - Output directory in which to save output files. Default is current - directory ('.'). + label : :obj:`str` or :obj:`None`, optional + Label for output directory. Default is None. Other Parameters ---------------- @@ -121,7 +120,10 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s', Notes ----- - This workflow writes out several files to ``out_dir``. + This workflow writes out several files, which are written out to a folder + named TED.[ref_label].[label] if ``label`` is provided and TED.[ref_label] + if not. ``ref_label`` is determined based on the name of the first ``data`` + file. Files are listed below: diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index 0a6b8981a..501b2d005 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -140,11 +140,11 @@ def _get_parser(): help='Perform PCA on wavelet-transformed data', action='store_true', default=False) - parser.add_argument('--out', - dest='out_dir', + parser.add_argument('--label', + dest='label', type=str, - help='Output directory.', - default='.') + help='Label for output directory.', + default=None) parser.add_argument('--seed', dest='fixed_seed', type=int, @@ -170,7 +170,7 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, kdaw=10., rdaw=1., conv=2.5e-5, ste=-1, combmode='t2s', dne=False, cost='logcosh', stabilize=False, filecsdata=False, wvpca=False, - out_dir='.', fixed_seed=42, debug=False, quiet=False): + label=None, fixed_seed=42, debug=False, quiet=False): """ Run the "canonical" TE-Dependent ANAlysis workflow. @@ -224,9 +224,8 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, wvpca : :obj:`bool`, optional Whether or not to perform PCA on wavelet-transformed data. Default is False. - out_dir : :obj:`str`, optional - Output directory in which to save output files. Default is current - directory ('.'). + label : :obj:`str` or :obj:`None`, optional + Label for output directory. Default is None. Other Parameters ---------------- @@ -291,6 +290,8 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, elif ctab is not None: raise IOError('Argument "ctab" must be an existing file.') + os.chdir(out_dir) + if mask is None: LGR.info('Computing adaptive mask') else: @@ -305,7 +306,7 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, # set a hard cap for the T2* map # anything that is 10x higher than the 99.5 %ile will be reset to 99.5 %ile - cap_t2s = stats.scoreatpercentile(t2s_limited.flatten(), 99.5, + cap_t2s = stats.scoreatpercentile(t2s.flatten(), 99.5, interpolation_method='lower') LGR.debug('Setting cap on T2* map at {:.5f}'.format(cap_t2s * 10)) t2s[t2s > cap_t2s * 10] = cap_t2s @@ -349,10 +350,10 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, else: LGR.info('Using supplied mixing matrix from ICA') mmix_orig = np.loadtxt(op.join(out_dir, 'meica_mix.1D')) - (seldict, comptable, - betas, mmix) = model.fitmodels_direct(catd, mmix_orig, - mask, t2s_limited, t2s_full, - tes, combmode, ref_img) + seldict, comptable, betas, mmix = model.fitmodels_direct(catd, mmix_orig, + mask, t2s, t2sG, + tes, combmode, + ref_img) if ctab is None: comptable = selection.selcomps(seldict, comptable, mmix, manacc, n_echos) From 08d4fb9c0554a60c602b5bf774afec71ca4df857 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 28 Nov 2018 18:20:56 -0500 Subject: [PATCH 19/21] Revert more. --- tedana/io.py | 16 ++++++---------- tedana/model/fit.py | 1 - tedana/workflows/t2smap.py | 2 +- 3 files changed, 7 insertions(+), 12 deletions(-) diff --git a/tedana/io.py b/tedana/io.py index fbdc7e88b..36181f309 100644 --- a/tedana/io.py +++ b/tedana/io.py @@ -273,7 +273,7 @@ def writefeats(data, mmix, mask, ref_img, suffix=''): def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, - acc, rej, midk, ign, ref_img): + acc, rej, midk, empty, ref_img): """ Denoises `ts` and saves all resulting files to disk @@ -300,7 +300,7 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, Indices of rejected (non-BOLD) components in `mmix` midk : :obj:`list` Indices of mid-K (questionable) components in `mmix` - ign : :obj:`list` + empty : :obj:`list` Indices of ignored components in `mmix` ref_img : :obj:`str` or img_like Reference image to dictate how outputs are saved to disk @@ -331,24 +331,20 @@ def writeresults(ts, mask, comptable, mmix, n_vols, fixed_seed, """ fout = filewrite(ts, 'ts_OC', ref_img) - LGR.info('Writing optimally combined time series: ' - '{}'.format(op.abspath(fout))) + LGR.info('Writing optimally combined time series: {}'.format(op.abspath(fout))) write_split_ts(ts, mmix, mask, acc, rej, midk, ref_img, suffix='OC') ts_B = model.get_coeffs(ts, mmix, mask) fout = filewrite(ts_B, 'betas_OC', ref_img) - LGR.info('Writing full ICA coefficient feature set: ' - '{}'.format(op.abspath(fout))) + LGR.info('Writing full ICA coefficient feature set: {}'.format(op.abspath(fout))) if len(acc) != 0: fout = filewrite(ts_B[:, acc], 'betas_hik_OC', ref_img) - LGR.info('Writing denoised ICA coefficient feature set: ' - '{}'.format(op.abspath(fout))) + LGR.info('Writing denoised ICA coefficient feature set: {}'.format(op.abspath(fout))) fout = writefeats(split_ts(ts, mmix, mask, acc)[0], mmix[:, acc], mask, ref_img, suffix='OC2') - LGR.info('Writing Z-normalized spatial component maps: ' - '{}'.format(op.abspath(fout))) + LGR.info('Writing Z-normalized spatial component maps: {}'.format(op.abspath(fout))) def writeresults_echoes(catd, mmix, mask, acc, rej, midk, ref_img): diff --git a/tedana/model/fit.py b/tedana/model/fit.py index 61b624db7..88302b005 100644 --- a/tedana/model/fit.py +++ b/tedana/model/fit.py @@ -2,7 +2,6 @@ Fit models. """ import logging -import os.path as op import numpy as np import pandas as pd diff --git a/tedana/workflows/t2smap.py b/tedana/workflows/t2smap.py index 156719945..85f901e24 100644 --- a/tedana/workflows/t2smap.py +++ b/tedana/workflows/t2smap.py @@ -108,7 +108,7 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s', combmode : {'t2s', 'ste'}, optional Combination scheme for TEs: 't2s' (Posse 1999, default), 'ste' (Poser). label : :obj:`str` or :obj:`None`, optional - Label for output directory. Default is None. + Label for output directory. Default is None. Other Parameters ---------------- From bb04f4ca0a48561933e57d77b82798e7f4646eb6 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 28 Nov 2018 18:23:50 -0500 Subject: [PATCH 20/21] Fix style problem. --- tedana/workflows/tedana.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py index 501b2d005..f8c16d9da 100644 --- a/tedana/workflows/tedana.py +++ b/tedana/workflows/tedana.py @@ -381,7 +381,7 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None, io.gscontrol_mmix(data_oc, mmix, mask, acc, ref_img, out_dir=out_dir) elif ws_denoise == 'godec': - decomposition.tedgodec(data_oc, mmix, mask, acc, empty, ref_img, + decomposition.tedgodec(data_oc, mmix, mask, acc, ign, ref_img, ranks=[2], wavelet=wvpca, thresh=10, norm_mode='vn', power=2, out_dir=out_dir) From 8b03fbb36892bfac19a9ceabbaf9efd6d39d4520 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 15 Feb 2019 09:37:53 -0500 Subject: [PATCH 21/21] Fix style issues. --- tedana/decomposition/__init__.py | 2 +- tedana/decomposition/godecomp.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tedana/decomposition/__init__.py b/tedana/decomposition/__init__.py index 3dd102614..ac77f5054 100644 --- a/tedana/decomposition/__init__.py +++ b/tedana/decomposition/__init__.py @@ -11,4 +11,4 @@ __all__ = [ 'tedpca', 'tedica', 'tedgodec', - ] +] diff --git a/tedana/decomposition/godecomp.py b/tedana/decomposition/godecomp.py index f239991fe..c159f8f61 100644 --- a/tedana/decomposition/godecomp.py +++ b/tedana/decomposition/godecomp.py @@ -52,7 +52,7 @@ def godec(data, thresh=.03, rank=2, power=1, tol=1e-3, max_iter=100, random_state = np.random.RandomState(random_seed) while True: Y2 = random_state.randn(n_vols, rank) - for i in range(power+1): + for i in range(power + 1): Y1 = np.dot(L, Y2) Y2 = np.dot(L.T, Y1) Q, R = qr(Y2) @@ -64,12 +64,12 @@ def godec(data, thresh=.03, rank=2, power=1, tol=1e-3, max_iter=100, err = np.linalg.norm(T.ravel(), 2) if err < tol: if verbose: - LGR.info('Successful convergence after %i iterations', itr+1) + LGR.info('Successful convergence after %i iterations', itr + 1) break elif itr >= max_iter: if verbose: LGR.warning('Model failed to converge after %i iterations', - itr+1) + itr + 1) break L += T itr += 1 @@ -106,7 +106,7 @@ def _tedgodec(data, wavelet=False, rank=2, power=2, tol=1e-3, # GoDec if wavelet: data_wt, cal = dwtmat(data_norm) - L, S, G = godec(data_wt, thresh=data_wt.std()*thresh, rank=rank, + L, S, G = godec(data_wt, thresh=(data_wt.std() * thresh), rank=rank, power=power, tol=tol, max_iter=max_iter, random_seed=random_seed, verbose=verbose) L = idwtmat(L, cal)