[REF] Add out_dir argument to filewriting functions (#500)

* Add out_dir argument to gscontrol functions.

* Add out_dir argument to io functions.

* Remove unnecessary chdir call.

* Fix removed files.

* Fix outputs.

* Fix io.filewrite.

* Fix filewrite.

tedana/gscontrol.py

@@ -2,6 +2,7 @@
 Global signal control methods
 """
 import logging
+import os.path as op
 
 import numpy as np
 from numpy.linalg import lstsq
@@ -15,7 +16,7 @@
 RefLGR = logging.getLogger('REFERENCES')
 
 
-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
 
@@ -35,6 +36,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`, optional
+        Output directory.
     dtrank : :obj:`int`, optional
         Specifies degree of Legendre polynomial basis function for estimating
         spatial global signal. Default: 4
@@ -75,23 +78,23 @@ 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()
-    io.filewrite(utils.unmask(sphis, Gmask), 'T1gs', ref_img)
+    io.filewrite(utils.unmask(sphis, Gmask), op.join(out_dir, '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('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
     sol = np.linalg.lstsq(np.atleast_2d(glbase), dat.T, rcond=None)[0]
     tsoc_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T,
                              np.atleast_2d(glbase.T[dtrank])) + Gmu[Gmask][:, np.newaxis]
 
-    io.filewrite(optcom, 'tsoc_orig', ref_img)
+    io.filewrite(optcom, op.join(out_dir, 'tsoc_orig'), ref_img)
     dm_optcom = utils.unmask(tsoc_nogs, Gmask)
-    io.filewrite(dm_optcom, 'tsoc_nogs', ref_img)
+    io.filewrite(dm_optcom, op.join(out_dir, 'tsoc_nogs'), ref_img)
 
     # Project glbase out of each echo
     dm_catd = catd.copy()  # don't overwrite catd
@@ -105,7 +108,7 @@ def gscontrol_raw(catd, optcom, n_echos, ref_img, dtrank=4):
     return dm_catd, dm_optcom
 
 
-def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img):
+def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img, out_dir='.'):
     """
     Perform global signal regression.
 
@@ -123,6 +126,8 @@ def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img):
         each metric. The index should be the component number.
     ref_img : :obj:`str` or img_like
         Reference image to dictate how outputs are saved to disk
+    out_dir : :obj:`str`, optional
+        Output directory.
 
     Notes
     -----
@@ -165,7 +170,7 @@ def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img):
     bold_ts = np.dot(cbetas[:, acc], mmix[:, acc].T)
     t1_map = bold_ts.min(axis=-1)
     t1_map -= t1_map.mean()
-    io.filewrite(utils.unmask(t1_map, mask), 'sphis_hik', ref_img)
+    io.filewrite(utils.unmask(t1_map, mask), op.join(out_dir, 'sphis_hik'), ref_img)
     t1_map = t1_map[:, np.newaxis]
 
     """
@@ -179,13 +184,14 @@ def gscontrol_mmix(optcom_ts, mmix, mask, comptable, 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
-    io.filewrite(utils.unmask(hik_ts, mask), 'hik_ts_OC_T1c.nii', ref_img)
+    io.filewrite(utils.unmask(hik_ts, mask), op.join(out_dir, 'hik_ts_OC_T1c'),
+                 ref_img)
 
     """
     Make denoised version of T1-corrected time series
     """
     medn_ts = optcom_mu + ((bold_noT1gs + resid) * optcom_std)
-    io.filewrite(utils.unmask(medn_ts, mask), 'dn_ts_OC_T1c.nii', ref_img)
+    io.filewrite(utils.unmask(medn_ts, mask), op.join(out_dir, 'dn_ts_OC_T1c'), ref_img)
 
     """
     Orthogonalize mixing matrix w.r.t. T1-GS
@@ -203,5 +209,5 @@ def gscontrol_mmix(optcom_ts, mmix, mask, comptable, ref_img):
     """
     cbetas_norm = lstsq(mmixnogs_norm.T, data_norm.T, rcond=None)[0].T
     io.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'), ref_img)
+    np.savetxt(op.join(out_dir, 'meica_mix_T1c.1D'), mmixnogs)

tedana/io.py

@@ -60,7 +60,7 @@ def split_ts(data, mmix, mask, comptable):
     return hikts, resid
 
 
-def write_split_ts(data, mmix, mask, comptable, ref_img, suffix=''):
+def write_split_ts(data, mmix, mask, comptable, ref_img, out_dir='.', suffix=''):
     """
     Splits `data` into denoised / noise / ignored time series and saves to disk
 
@@ -75,6 +75,8 @@ def write_split_ts(data, mmix, mask, comptable, ref_img, suffix=''):
         Boolean mask array
     ref_img : :obj:`str` or img_like
         Reference image to dictate how outputs are saved to disk
+    out_dir : :obj:`str`, optional
+        Output directory.
     suffix : :obj:`str`, optional
         Appended to name of saved files (before extension). Default: ''
 
@@ -116,22 +118,21 @@ def write_split_ts(data, mmix, mask, comptable, ref_img, suffix=''):
 
     if len(acc) != 0:
         fout = filewrite(utils.unmask(hikts, mask),
-                         'hik_ts_{0}'.format(suffix), ref_img)
+                         op.join(out_dir, 'hik_ts_{0}'.format(suffix)), ref_img)
         LGR.info('Writing high-Kappa time series: {}'.format(op.abspath(fout)))
 
     if len(rej) != 0:
         fout = filewrite(utils.unmask(lowkts, mask),
-                         'lowk_ts_{0}'.format(suffix), ref_img)
+                         op.join(out_dir, 'lowk_ts_{0}'.format(suffix)), ref_img)
         LGR.info('Writing low-Kappa time series: {}'.format(op.abspath(fout)))
 
     fout = filewrite(utils.unmask(dnts, mask),
-                     'dn_ts_{0}'.format(suffix), ref_img)
+                     op.join(out_dir, '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=''):
+def writefeats(data, mmix, mask, ref_img, out_dir='.', suffix=''):
     """
     Converts `data` to component space with `mmix` and saves to disk
 
@@ -146,6 +147,8 @@ def writefeats(data, mmix, mask, ref_img, suffix=''):
         Boolean mask array
     ref_img : :obj:`str` or img_like
         Reference image to dictate how outputs are saved to disk
+    out_dir : :obj:`str`, optional
+        Output directory.
     suffix : :obj:`str`, optional
         Appended to name of saved files (before extension). Default: ''
 
@@ -167,12 +170,11 @@ def writefeats(data, mmix, mask, ref_img, suffix=''):
 
     # write feature versions of components
     feats = utils.unmask(computefeats2(data, mmix, mask), mask)
-    fname = filewrite(feats, 'feats_{0}'.format(suffix), ref_img)
-
+    fname = filewrite(feats, op.join(out_dir, 'feats_{0}'.format(suffix)), ref_img)
     return fname
 
 
-def writeresults(ts, mask, comptable, mmix, n_vols, ref_img):
+def writeresults(ts, mask, comptable, mmix, n_vols, ref_img, out_dir='.'):
     """
     Denoises `ts` and saves all resulting files to disk
 
@@ -193,6 +195,8 @@ def writeresults(ts, mask, comptable, mmix, n_vols, ref_img):
         Number of volumes in original time series
     ref_img : :obj:`str` or img_like
         Reference image to dictate how outputs are saved to disk
+    out_dir : :obj:`str`, optional
+        Output directory.
 
     Notes
     -----
@@ -201,7 +205,6 @@ def writeresults(ts, mask, comptable, mmix, n_vols, ref_img):
     ======================    =================================================
     Filename                  Content
     ======================    =================================================
-    ts_OC.nii                 Optimally combined 4D time series.
     hik_ts_OC.nii             High-Kappa time series. Generated by
                               :py:func:`tedana.utils.io.write_split_ts`.
     midk_ts_OC.nii            Mid-Kappa time series. Generated by
@@ -214,28 +217,30 @@ def writeresults(ts, mask, comptable, mmix, n_vols, ref_img):
     betas_hik_OC.nii          Denoised ICA coefficient feature set.
     feats_OC2.nii             Z-normalized spatial component maps. Generated
                               by :py:func:`tedana.utils.io.writefeats`.
+    ts_OC.nii                 Optimally combined 4D time series.
     ======================    =================================================
     """
     acc = comptable[comptable.classification == 'accepted'].index.values
 
-    fout = filewrite(ts, 'ts_OC', ref_img)
+    fout = filewrite(ts, op.join(out_dir, 'ts_OC'), ref_img)
     LGR.info('Writing optimally combined time series: {}'.format(op.abspath(fout)))
 
-    write_split_ts(ts, mmix, mask, comptable, ref_img, suffix='OC')
+    write_split_ts(ts, mmix, mask, comptable, ref_img, out_dir=out_dir, suffix='OC')
 
     ts_B = get_coeffs(ts, mmix, mask)
-    fout = filewrite(ts_B, 'betas_OC', ref_img)
+    fout = filewrite(ts_B, op.join(out_dir, '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], 'betas_hik_OC', ref_img)
+        fout = filewrite(ts_B[:, acc], op.join(out_dir, 'betas_hik_OC'), ref_img)
         LGR.info('Writing denoised ICA coefficient feature set: {}'.format(op.abspath(fout)))
         fout = writefeats(split_ts(ts, mmix, mask, comptable)[0],
-                          mmix[:, acc], mask, ref_img, suffix='OC2')
+                          mmix[:, acc], mask, ref_img, out_dir=out_dir,
+                          suffix='OC2')
         LGR.info('Writing Z-normalized spatial component maps: {}'.format(op.abspath(fout)))
 
 
-def writeresults_echoes(catd, mmix, mask, comptable, ref_img):
+def writeresults_echoes(catd, mmix, mask, comptable, ref_img, out_dir='.'):
     """
     Saves individually denoised echos to disk
 
@@ -253,6 +258,8 @@ def writeresults_echoes(catd, mmix, mask, comptable, ref_img):
         each metric. The index should be the component number.
     ref_img : :obj:`str` or img_like
         Reference image to dictate how outputs are saved to disk
+    out_dir : :obj:`str`, optional
+        Output directory.
 
     Notes
     -----
@@ -279,7 +286,7 @@ def writeresults_echoes(catd, mmix, mask, comptable, 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, comptable, ref_img,
-                       suffix='e%i' % (i_echo + 1))
+                       out_dir=out_dir, suffix='e%i' % (i_echo + 1))
 
 
 def new_nii_like(ref_img, data, affine=None, copy_header=True):
@@ -351,7 +358,10 @@ def filewrite(data, filename, ref_img, gzip=True, copy_header=True):
 
     # FIXME: we only handle writing to nifti right now
     # get root of desired output file and save as nifti image
-    root, ext, add = splitext_addext(filename)
+    root = op.dirname(filename)
+    base = op.basename(filename)
+    base, ext, add = splitext_addext(base)
+    root = op.join(root, base)
     name = '{}.{}'.format(root, 'nii.gz' if gzip else 'nii')
     out.to_filename(name)
 

tedana/workflows/tedana.py

@@ -462,8 +462,6 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
     LGR.debug('Retaining {}/{} samples'.format(mask.sum(), n_samp))
     io.filewrite(masksum, op.join(out_dir, 'adaptive_mask.nii'), ref_img)
 
-    os.chdir(out_dir)
-
     if t2smap is None:
         LGR.info('Computing T2* map')
         t2s_limited, s0_limited, t2s_full, s0_full = decay.fit_decay(
@@ -487,7 +485,8 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
 
     # regress out global signal unless explicitly not desired
     if 'gsr' in gscontrol:
-        catd, data_oc = gsc.gscontrol_raw(catd, data_oc, n_echos, ref_img)
+        catd, data_oc = gsc.gscontrol_raw(catd, data_oc, n_echos, ref_img,
+                                          out_dir=out_dir)
 
     if mixm is None:
         # Identify and remove thermal noise from data
@@ -516,7 +515,7 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
         comp_names = [io.add_decomp_prefix(comp, prefix='ica', max_value=comptable.index.max())
                       for comp in comptable.index.values]
         mixing_df = pd.DataFrame(data=mmix, columns=comp_names)
-        mixing_df.to_csv('ica_mixing.tsv', sep='\t', index=False)
+        mixing_df.to_csv(op.join(out_dir, 'ica_mixing.tsv'), sep='\t', index=False)
         betas_oc = utils.unmask(computefeats2(data_oc, mmix, mask), mask)
         io.filewrite(betas_oc,
                      op.join(out_dir, 'ica_components.nii.gz'),
@@ -575,7 +574,7 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
         comp_names = [io.add_decomp_prefix(comp, prefix='ica', max_value=comptable.index.max())
                       for comp in comptable.index.values]
         mixing_df = pd.DataFrame(data=mmix, columns=comp_names)
-        mixing_df.to_csv('ica_orth_mixing.tsv', sep='\t', index=False)
+        mixing_df.to_csv(op.join(out_dir, 'ica_orth_mixing.tsv'), sep='\t', index=False)
         RepLGR.info("Rejected components' time series were then "
                     "orthogonalized with respect to accepted components' time "
                     "series.")
@@ -585,13 +584,14 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
                     comptable=comptable,
                     mmix=mmix,
                     n_vols=n_vols,
-                    ref_img=ref_img)
+                    ref_img=ref_img,
+                    out_dir=out_dir)
 
     if 't1c' in gscontrol:
-        gsc.gscontrol_mmix(data_oc, mmix, mask, comptable, ref_img)
+        gsc.gscontrol_mmix(data_oc, mmix, mask, comptable, ref_img, out_dir=out_dir)
 
     if verbose:
-        io.writeresults_echoes(catd, mmix, mask, comptable, ref_img)
+        io.writeresults_echoes(catd, mmix, mask, comptable, ref_img, out_dir=out_dir)
 
     if not no_png:
         LGR.info('Making figures folder with static component maps and '