From 86a73e1db0c2a9db4c07e4371e7c86da083c16e4 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 10 Apr 2024 13:49:58 -0400
Subject: [PATCH 01/12] Create denoise_echoes.py

---
 tedana/workflows/denoise_echoes.py | 414 +++++++++++++++++++++++++++++
 1 file changed, 414 insertions(+)
 create mode 100644 tedana/workflows/denoise_echoes.py

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
new file mode 100644
index 000000000..ebc27d1f0
--- /dev/null
+++ b/tedana/workflows/denoise_echoes.py
@@ -0,0 +1,414 @@
+"""Estimate T2 and S0, and optimally combine data across TEs."""
+
+import argparse
+import logging
+import os
+import os.path as op
+import sys
+
+import numpy as np
+from scipy import stats
+from threadpoolctl import threadpool_limits
+
+from tedana import __version__, combine, decay, io, utils
+from tedana.workflows.parser_utils import is_valid_file
+
+LGR = logging.getLogger("GENERAL")
+RepLGR = logging.getLogger("REPORT")
+
+
+def _get_parser():
+    """Parse command line inputs for denoise_echoes.
+
+    Returns
+    -------
+    parser.parse_args() : argparse dict
+    """
+    parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    # Argument parser follow template provided by RalphyZ
+    # https://stackoverflow.com/a/43456577
+    optional = parser._action_groups.pop()
+    required = parser.add_argument_group("Required Arguments")
+    required.add_argument(
+        "-d",
+        dest="data",
+        nargs="+",
+        metavar="FILE",
+        type=lambda x: is_valid_file(parser, x),
+        help=(
+            "Multi-echo dataset for analysis. May be a "
+            "single file with spatially concatenated data "
+            "or a set of echo-specific files, in the same "
+            "order as the TEs are listed in the -e "
+            "argument."
+        ),
+        required=True,
+    )
+    required.add_argument(
+        "-e",
+        dest="tes",
+        nargs="+",
+        metavar="TE",
+        type=float,
+        help="Echo times (in ms). E.g., 15.0 39.0 63.0",
+        required=True,
+    )
+    required.add_argument(
+        "--confounds",
+        dest="confounds",
+        nargs="+",
+        metavar="FILE",
+        type=lambda x: is_valid_file(parser, x),
+        help="Files defining confounds to regress from the echo-wise data.",
+    )
+    optional.add_argument(
+        "--out-dir",
+        dest="out_dir",
+        type=str,
+        metavar="PATH",
+        help="Output directory.",
+        default=".",
+    )
+    optional.add_argument(
+        "--mask",
+        dest="mask",
+        metavar="FILE",
+        type=lambda x: is_valid_file(parser, x),
+        help=(
+            "Binary mask of voxels to include in TE "
+            "Dependent ANAlysis. Must be in the same "
+            "space as `data`."
+        ),
+        default=None,
+    )
+    optional.add_argument(
+        "--prefix", dest="prefix", type=str, help="Prefix for filenames generated.", default=""
+    )
+    optional.add_argument(
+        "--convention",
+        dest="convention",
+        action="store",
+        choices=["orig", "bids"],
+        help=("Filenaming convention. bids will use the latest BIDS derivatives version."),
+        default="bids",
+    )
+    optional.add_argument(
+        "--masktype",
+        dest="masktype",
+        required=False,
+        action="store",
+        nargs="+",
+        help="Method(s) by which to define the adaptive mask.",
+        choices=["dropout", "decay", "none"],
+        default=["dropout"],
+    )
+    optional.add_argument(
+        "--fittype",
+        dest="fittype",
+        action="store",
+        choices=["loglin", "curvefit"],
+        help=(
+            "Desired T2*/S0 fitting method. "
+            '"loglin" means that a linear model is fit '
+            "to the log of the data. "
+            '"curvefit" means that a more computationally '
+            "demanding monoexponential model is fit "
+            "to the raw data. "
+        ),
+        default="loglin",
+    )
+    optional.add_argument(
+        "--fitmode",
+        dest="fitmode",
+        action="store",
+        choices=["all", "ts"],
+        help=(
+            "Monoexponential model fitting scheme. "
+            '"all" means that the model is fit, per voxel, '
+            "across all timepoints. "
+            '"ts" means that the model is fit, per voxel '
+            "and per timepoint."
+        ),
+        default="all",
+    )
+    optional.add_argument(
+        "--combmode",
+        dest="combmode",
+        action="store",
+        choices=["t2s", "paid"],
+        help=("Combination scheme for TEs: t2s (Posse 1999), paid (Poser)"),
+        default="t2s",
+    )
+    optional.add_argument(
+        "--n-threads",
+        dest="n_threads",
+        type=int,
+        action="store",
+        help=(
+            "Number of threads to use. Used by "
+            "threadpoolctl to set the parameter outside "
+            "of the workflow function. Higher numbers of "
+            "threads tend to slow down performance on "
+            "typical datasets."
+        ),
+        default=1,
+    )
+    optional.add_argument(
+        "--debug", dest="debug", help=argparse.SUPPRESS, action="store_true", default=False
+    )
+    optional.add_argument(
+        "--quiet", dest="quiet", help=argparse.SUPPRESS, action="store_true", default=False
+    )
+    parser._action_groups.append(optional)
+    return parser
+
+
+def t2smap_workflow(
+    data,
+    tes,
+    out_dir=".",
+    confounds=None,
+    mask=None,
+    prefix="",
+    convention="bids",
+    masktype=["dropout"],
+    fittype="loglin",
+    fitmode="all",
+    combmode="t2s",
+    debug=False,
+    quiet=False,
+    t2smap_command=None,
+):
+    """
+    Estimate T2 and S0, and optimally combine data across TEs.
+
+    Please remember to cite :footcite:t:`dupre2021te`.
+
+    Parameters
+    ----------
+    data : :obj:`str` or :obj:`list` of :obj:`str`
+        Either a single z-concatenated file (single-entry list or str) or a
+        list of echo-specific files, in ascending order.
+    tes : :obj:`list`
+        List of echo times associated with data in milliseconds.
+    out_dir : :obj:`str`, optional
+        Output directory.
+    confounds : :obj:`list` of :obj:`str`
+        Files defining confounds to regress from the echo-wise data.
+    mask : :obj:`str`, optional
+        Binary mask of voxels to include in TE Dependent ANAlysis. Must be spatially
+        aligned with `data`.
+    masktype : :obj:`list` with 'dropout' and/or 'decay' or None, optional
+        Method(s) by which to define the adaptive mask. Default is ["dropout"].
+    fittype : {'loglin', 'curvefit'}, optional
+        Monoexponential fitting method.
+        'loglin' means to use the the default linear fit to the log of
+        the data.
+        'curvefit' means to use a monoexponential fit to the raw data,
+        which is slightly slower but may be more accurate.
+    fitmode : {'all', 'ts'}, optional
+        Monoexponential model fitting scheme.
+        'all' means that the model is fit, per voxel, across all timepoints.
+        'ts' means that the model is fit, per voxel and per timepoint.
+        Default is 'all'.
+    combmode : {'t2s', 'paid'}, optional
+        Combination scheme for TEs: 't2s' (Posse 1999, default), 'paid' (Poser).
+    t2smap_command : :obj:`str`, optional
+        The command used to run t2smap. Default is None.
+
+    Other Parameters
+    ----------------
+    debug : :obj:`bool`, optional
+        Whether to run in debugging mode or not. Default is False.
+    quiet : :obj:`bool`, optional
+        If True, suppress logging/printing of messages. Default is False.
+
+    Notes
+    -----
+    This workflow writes out several files, which are described below:
+
+    ============================= =================================================
+    Filename                      Content
+    ============================= =================================================
+    T2starmap.nii.gz              Estimated T2* 3D map or 4D timeseries.
+                                  Will be a 3D map if ``fitmode`` is 'all' and a
+                                  4D timeseries if it is 'ts'.
+    S0map.nii.gz                  S0 3D map or 4D timeseries.
+    desc-limited_T2starmap.nii.gz Limited T2* map/timeseries. The difference between
+                                  the limited and full maps is that, for voxels
+                                  affected by dropout where only one echo contains
+                                  good data, the full map uses the T2* estimate
+                                  from the first two echos, while the limited map
+                                  will have a NaN.
+    desc-limited_S0map.nii.gz     Limited S0 map/timeseries. The difference between
+                                  the limited and full maps is that, for voxels
+                                  affected by dropout where only one echo contains
+                                  good data, the full map uses the S0 estimate
+                                  from the first two echos, while the limited map
+                                  will have a NaN.
+    desc-optcom_bold.nii.gz       Optimally combined timeseries with the confounds
+                                  regressed out.
+    echo-<n>_bold.nii.gz          Echo-specific timeseries with the confounds
+                                  regressed out.
+    ============================= =================================================
+
+    References
+    ----------
+    .. footbibliography::
+    """
+    out_dir = op.abspath(out_dir)
+    if not op.isdir(out_dir):
+        os.mkdir(out_dir)
+
+    utils.setup_loggers(quiet=quiet, debug=debug)
+
+    LGR.info(f"Using output directory: {out_dir}")
+
+    # Save command into sh file, if the command-line interface was used
+    if t2smap_command is not None:
+        command_file = open(os.path.join(out_dir, "t2smap_call.sh"), "w")
+        command_file.write(t2smap_command)
+        command_file.close()
+    else:
+        # Get variables passed to function if the tedana command is None
+        variables = ", ".join(f"{name}={value}" for name, value in locals().items())
+        # From variables, remove everything after ", tedana_command"
+        variables = variables.split(", t2smap_command")[0]
+        t2smap_command = f"t2smap_workflow({variables})"
+
+    # Save system info to json
+    info_dict = utils.get_system_version_info()
+    info_dict["Command"] = t2smap_command
+
+    # ensure tes are in appropriate format
+    tes = [float(te) for te in tes]
+    n_echos = len(tes)
+
+    # coerce data to samples x echos x time array
+    if isinstance(data, str):
+        data = [data]
+
+    LGR.info(f"Loading input data: {[f for f in data]}")
+    catd, ref_img = io.load_data(data, n_echos=n_echos)
+    io_generator = io.OutputGenerator(
+        ref_img,
+        convention=convention,
+        out_dir=out_dir,
+        prefix=prefix,
+        config="auto",
+        make_figures=False,
+    )
+    n_samp, n_echos, n_vols = catd.shape
+    LGR.debug(f"Resulting data shape: {catd.shape}")
+
+    if mask is None:
+        LGR.info("Computing adaptive mask")
+    else:
+        LGR.info("Using user-defined mask")
+    mask, masksum = utils.make_adaptive_mask(
+        catd,
+        mask=mask,
+        threshold=1,
+        methods=masktype,
+    )
+
+    LGR.info("Computing adaptive T2* map")
+    if fitmode == "all":
+        (t2s_limited, s0_limited, t2s_full, s0_full) = decay.fit_decay(
+            catd, tes, mask, masksum, fittype
+        )
+    else:
+        (t2s_limited, s0_limited, t2s_full, s0_full) = decay.fit_decay_ts(
+            catd, tes, mask, masksum, fittype
+        )
+
+    # set a hard cap for the T2* map/timeseries
+    # anything that is 10x higher than the 99.5 %ile will be reset to 99.5 %ile
+    cap_t2s = stats.scoreatpercentile(t2s_full.flatten(), 99.5, interpolation_method="lower")
+    cap_t2s_sec = utils.millisec2sec(cap_t2s * 10.0)
+    LGR.debug(f"Setting cap on T2* map at {cap_t2s_sec:.5f}s")
+    t2s_full[t2s_full > cap_t2s * 10] = cap_t2s
+
+    LGR.info("Computing optimal combination")
+    # optimally combine data
+    data_oc = combine.make_optcom(catd, tes, masksum, t2s=t2s_full, combmode=combmode)
+
+    # clean up numerical errors
+    for arr in (data_oc, s0_full, t2s_full):
+        np.nan_to_num(arr, copy=False)
+
+    if "gsr" in gscontrol:
+        # regress out global signal
+        catd, data_oc = gsc.gscontrol_raw(catd, data_oc, n_echos, io_generator)
+
+    s0_full[s0_full < 0] = 0
+    t2s_full[t2s_full < 0] = 0
+
+    io_generator.save_file(
+        utils.millisec2sec(t2s_full),
+        "t2star img",
+    )
+    io_generator.save_file(s0_full, "s0 img")
+    io_generator.save_file(
+        utils.millisec2sec(t2s_limited),
+        "limited t2star img",
+    )
+    io_generator.save_file(
+        s0_limited,
+        "limited s0 img",
+    )
+    io_generator.save_file(data_oc, "combined img")
+
+    # Write out BIDS-compatible description file
+    derivative_metadata = {
+        "Name": "t2smap Outputs",
+        "BIDSVersion": "1.5.0",
+        "DatasetType": "derivative",
+        "GeneratedBy": [
+            {
+                "Name": "t2smap",
+                "Version": __version__,
+                "Description": (
+                    "A pipeline estimating T2* from multi-echo fMRI data and "
+                    "combining data across echoes."
+                ),
+                "CodeURL": "https://github.com/ME-ICA/tedana",
+                "Node": {
+                    "Name": info_dict["Node"],
+                    "System": info_dict["System"],
+                    "Machine": info_dict["Machine"],
+                    "Processor": info_dict["Processor"],
+                    "Release": info_dict["Release"],
+                    "Version": info_dict["Version"],
+                },
+                "Python": info_dict["Python"],
+                "Python_Libraries": info_dict["Python_Libraries"],
+                "Command": info_dict["Command"],
+            }
+        ],
+    }
+
+    io_generator.save_file(derivative_metadata, "data description json")
+    io_generator.save_self()
+
+    LGR.info("Workflow completed")
+    utils.teardown_loggers()
+
+
+def _main(argv=None):
+    """Run the t2smap workflow."""
+    if argv:
+        # relevant for tests when CLI called with t2smap_cli._main(args)
+        t2smap_command = "t2smap " + " ".join(argv)
+    else:
+        t2smap_command = "t2smap " + " ".join(sys.argv[1:])
+    options = _get_parser().parse_args(argv)
+    kwargs = vars(options)
+    n_threads = kwargs.pop("n_threads")
+    n_threads = None if n_threads == -1 else n_threads
+    with threadpool_limits(limits=n_threads, user_api=None):
+        t2smap_workflow(**kwargs, t2smap_command=t2smap_command)
+
+
+if __name__ == "__main__":
+    _main()

From ee0d37f8401d4e99f44ef51366e16a3c097ec91f Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Sun, 21 Apr 2024 13:15:24 -0400
Subject: [PATCH 02/12] Keep working.

---
 tedana/workflows/denoise_echoes.py | 108 +++++++++++++++++++++++++----
 1 file changed, 94 insertions(+), 14 deletions(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index ebc27d1f0..8c026f849 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -7,6 +7,7 @@
 import sys
 
 import numpy as np
+import pandas as pd
 from scipy import stats
 from threadpoolctl import threadpool_limits
 
@@ -163,7 +164,7 @@ def _get_parser():
     return parser
 
 
-def t2smap_workflow(
+def denoise_echoes_workflow(
     data,
     tes,
     out_dir=".",
@@ -289,7 +290,7 @@ def t2smap_workflow(
         data = [data]
 
     LGR.info(f"Loading input data: {[f for f in data]}")
-    catd, ref_img = io.load_data(data, n_echos=n_echos)
+    data_cat, ref_img = io.load_data(data, n_echos=n_echos)
     io_generator = io.OutputGenerator(
         ref_img,
         convention=convention,
@@ -298,15 +299,15 @@ def t2smap_workflow(
         config="auto",
         make_figures=False,
     )
-    n_samp, n_echos, n_vols = catd.shape
-    LGR.debug(f"Resulting data shape: {catd.shape}")
+    n_samp, n_echos, n_vols = data_cat.shape
+    LGR.debug(f"Resulting data shape: {data_cat.shape}")
 
     if mask is None:
         LGR.info("Computing adaptive mask")
     else:
         LGR.info("Using user-defined mask")
-    mask, masksum = utils.make_adaptive_mask(
-        catd,
+    mask, adaptive_mask = utils.make_adaptive_mask(
+        data_cat,
         mask=mask,
         threshold=1,
         methods=masktype,
@@ -315,11 +316,11 @@ def t2smap_workflow(
     LGR.info("Computing adaptive T2* map")
     if fitmode == "all":
         (t2s_limited, s0_limited, t2s_full, s0_full) = decay.fit_decay(
-            catd, tes, mask, masksum, fittype
+            data_cat, tes, mask, adaptive_mask, fittype
         )
     else:
         (t2s_limited, s0_limited, t2s_full, s0_full) = decay.fit_decay_ts(
-            catd, tes, mask, masksum, fittype
+            data_cat, tes, mask, adaptive_mask, fittype
         )
 
     # set a hard cap for the T2* map/timeseries
@@ -331,16 +332,12 @@ def t2smap_workflow(
 
     LGR.info("Computing optimal combination")
     # optimally combine data
-    data_oc = combine.make_optcom(catd, tes, masksum, t2s=t2s_full, combmode=combmode)
+    data_oc = combine.make_optcom(data_cat, tes, adaptive_mask, t2s=t2s_full, combmode=combmode)
 
     # clean up numerical errors
     for arr in (data_oc, s0_full, t2s_full):
         np.nan_to_num(arr, copy=False)
 
-    if "gsr" in gscontrol:
-        # regress out global signal
-        catd, data_oc = gsc.gscontrol_raw(catd, data_oc, n_echos, io_generator)
-
     s0_full[s0_full < 0] = 0
     t2s_full[t2s_full < 0] = 0
 
@@ -407,8 +404,91 @@ def _main(argv=None):
     n_threads = kwargs.pop("n_threads")
     n_threads = None if n_threads == -1 else n_threads
     with threadpool_limits(limits=n_threads, user_api=None):
-        t2smap_workflow(**kwargs, t2smap_command=t2smap_command)
+        denoise_echoes_workflow(**kwargs, t2smap_command=t2smap_command)
 
 
 if __name__ == "__main__":
     _main()
+
+
+def denoise_echoes(data_cat, data_oc, mask, confounds):
+    LGR.info("Applying amplitude-based T1 equilibration correction")
+    RepLGR.info(
+        "Global signal regression was applied to the multi-echo "
+        "and optimally combined datasets."
+    )
+    n_samples, n_echos, n_volumes = data_cat.shape
+    if data_cat.shape[0] != data_oc.shape[0]:
+        raise ValueError(
+            f"First dimensions of data_cat ({data_cat.shape[0]}) and data_oc ({data_oc.shape[0]}) "
+            "do not match"
+        )
+    elif data_cat.shape[2] != data_oc.shape[1]:
+        raise ValueError(
+            f"Third dimension of data_cat ({data_cat.shape[2]}) does not match second dimension "
+            f"of data_oc ({data_oc.shape[1]})"
+        )
+
+    confound_arrs = []
+    confound_types = []
+    for confound in confounds:
+        if confound.endswith(".tsv"):
+            # One or more time series
+            confound_df = pd.read_table(confound)
+            confound_arr = confound_df.to_numpy()
+            assert confound_arr.shape[0] == n_volumes
+            confound_types.append("temporal")
+        elif confound.endswith(".nii.gz"):
+            confound_arr = utils.reshape_niimg(confound)
+            if confound_arr.ndim == 1:
+                # Spatial map that must be regressed to produce a single time series
+                assert confound_arr.shape[0] == n_samples
+                confound_types.append("spatial")
+
+                # Mask out out-of-brain voxels
+                confound_arr = confound_arr[mask]
+
+                # Find the time course for the spatial map
+                confound_timeseries = np.linalg.lstsq(
+                    np.atleast_2d(confound_arr).T,
+                    data_oc,
+                    rcond=None,
+                )[0]
+                confound_timeseries = stats.zscore(confound_timeseries, axis=None)
+
+            elif confound_arr.ndim == 2:
+                # Voxel-wise regressors
+                assert confound_arr.shape[0] == n_samples
+                assert confound_arr.shape[1] == n_volumes
+                # Mask out out-of-brain voxels
+                confound_arr = confound_arr[mask, ...]
+                confound_types.append("voxelwise")
+
+            else:
+                raise ValueError(f"Unknown shape for confound array: {confound_arr.shape}")
+        else:
+            raise ValueError(f"Unknown file type: {confound}")
+
+        confound_arrs.append(confound_arr)
+
+    # 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]))
+        + temporal_mean[temporal_mean_mask][:, np.newaxis]
+    )
+
+    io_generator.save_file(data_oc, "has gs combined img")
+    dm_optcom = utils.unmask(tsoc_nogs, temporal_mean_mask)
+    io_generator.save_file(dm_optcom, "removed gs combined img")
+
+    # Project glbase out of each echo
+    dm_catd = data_cat.copy()  # don't overwrite data_cat
+    for echo in range(n_echos):
+        dat = dm_catd[:, echo, :][temporal_mean_mask]
+        sol = np.linalg.lstsq(np.atleast_2d(glbase), dat.T, rcond=None)[0]
+        e_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T, np.atleast_2d(glbase.T[dtrank]))
+        dm_catd[:, echo, :] = utils.unmask(e_nogs, temporal_mean_mask)
+
+    return dm_catd, dm_optcom

From d182ba9899d17023b3568ee3f1a1b1fdc303483c Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Mon, 22 Apr 2024 10:27:37 -0400
Subject: [PATCH 03/12] Update denoise_echoes.py

---
 tedana/workflows/denoise_echoes.py | 128 ++++++++++++++++++++---------
 1 file changed, 90 insertions(+), 38 deletions(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index 8c026f849..10fbdbf93 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -10,6 +10,7 @@
 import pandas as pd
 from scipy import stats
 from threadpoolctl import threadpool_limits
+from tqdm import trange
 
 from tedana import __version__, combine, decay, io, utils
 from tedana.workflows.parser_utils import is_valid_file
@@ -411,12 +412,37 @@ def _main(argv=None):
     _main()
 
 
-def denoise_echoes(data_cat, data_oc, mask, confounds):
-    LGR.info("Applying amplitude-based T1 equilibration correction")
-    RepLGR.info(
-        "Global signal regression was applied to the multi-echo "
-        "and optimally combined datasets."
-    )
+def denoise_echoes(
+    data_cat: np.ndarray,
+    data_oc: np.ndarray,
+    mask: np.ndarray,
+    confounds: dict[str, str],
+) -> tuple[np.ndarray, np.ndarray]:
+    """Denoise echoes using external regressors.
+
+    TODO: Calculate confound-wise echo-dependence metrics.
+
+    Parameters
+    ----------
+    data_cat : :obj:`numpy.ndarray` of shape (n_samples, n_echos, n_volumes)
+        Concatenated data across echoes.
+    data_oc : :obj:`numpy.ndarray` of shape (n_samples, n_volumes)
+        Optimally combined data across echoes.
+    mask : :obj:`numpy.ndarray` of shape (n_samples,)
+        Binary mask of voxels to include in TE Dependent ANAlysis.
+    confounds : :obj:`dict`
+        Files defining confounds to regress from the echo-wise data.
+        Keys indicate the names to use for the confounds.
+
+    Returns
+    -------
+    data_cat_denoised : :obj:`numpy.ndarray` of shape (n_samples, n_echos, n_volumes)
+        Denoised concatenated data across echoes.
+    data_optcom_denoised : :obj:`numpy.ndarray` of shape (n_samples, n_volumes)
+        Denoised optimally combined data.
+    """
+    LGR.info("Applying a priori confound regression")
+    RepLGR.info("Confounds were regressed out of the multi-echo and optimally combined datasets.")
     n_samples, n_echos, n_volumes = data_cat.shape
     if data_cat.shape[0] != data_oc.shape[0]:
         raise ValueError(
@@ -429,21 +455,20 @@ def denoise_echoes(data_cat, data_oc, mask, confounds):
             f"of data_oc ({data_oc.shape[1]})"
         )
 
-    confound_arrs = []
-    confound_types = []
-    for confound in confounds:
-        if confound.endswith(".tsv"):
+    confound_dfs = []
+    voxelwise_confounds = {}
+    for confound_name, confound_file in confounds.items():
+        if confound_file.endswith(".tsv"):
             # One or more time series
-            confound_df = pd.read_table(confound)
-            confound_arr = confound_df.to_numpy()
-            assert confound_arr.shape[0] == n_volumes
-            confound_types.append("temporal")
-        elif confound.endswith(".nii.gz"):
-            confound_arr = utils.reshape_niimg(confound)
+            confound_df = pd.read_table(confound_file)
+            confound_df = confound_df.add_suffix(f"_{confound_name}")
+            assert confound_df.shape[0] == n_volumes
+            confound_dfs.append(confound_df)
+        elif confound_file.endswith(".nii.gz"):
+            confound_arr = utils.reshape_niimg(confound_file)
             if confound_arr.ndim == 1:
                 # Spatial map that must be regressed to produce a single time series
                 assert confound_arr.shape[0] == n_samples
-                confound_types.append("spatial")
 
                 # Mask out out-of-brain voxels
                 confound_arr = confound_arr[mask]
@@ -455,6 +480,9 @@ def denoise_echoes(data_cat, data_oc, mask, confounds):
                     rcond=None,
                 )[0]
                 confound_timeseries = stats.zscore(confound_timeseries, axis=None)
+                confound_df = pd.DataFrame(confound_timeseries.T, columns=[confound_name])
+                assert confound_df.shape[0] == n_volumes
+                confound_dfs.append(confound_df)
 
             elif confound_arr.ndim == 2:
                 # Voxel-wise regressors
@@ -462,33 +490,57 @@ def denoise_echoes(data_cat, data_oc, mask, confounds):
                 assert confound_arr.shape[1] == n_volumes
                 # Mask out out-of-brain voxels
                 confound_arr = confound_arr[mask, ...]
-                confound_types.append("voxelwise")
+                voxelwise_confounds[confound_name] = confound_arr
 
             else:
                 raise ValueError(f"Unknown shape for confound array: {confound_arr.shape}")
         else:
-            raise ValueError(f"Unknown file type: {confound}")
+            raise ValueError(f"Unknown file type: {confound_file}")
 
-        confound_arrs.append(confound_arr)
-
-    # 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]))
-        + temporal_mean[temporal_mean_mask][:, np.newaxis]
-    )
+    if confound_dfs:
+        confounds_df = pd.concat(confound_dfs, axis=1)
+        confounds_df = confounds_df.fillna(0)
+    else:
+        confounds_df = pd.DataFrame(index=range(n_volumes))
+
+    # Project confounds out of optimally combined data
+    temporal_mean = data_oc.mean(axis=-1)  # temporal mean
+    data_optcom_denoised = data_oc[mask] - temporal_mean[mask, np.newaxis]
+    if voxelwise_confounds:
+        for i_voxel in trange(data_optcom_denoised.shape[0], desc="Denoise optimally combined"):
+            design_matrix = confounds_df.copy()
+            for confound_name, confound_arr in voxelwise_confounds.items():
+                design_matrix[confound_name] = confound_arr[i_voxel, :]
+
+            betas = np.linalg.lstsq(design_matrix.values, data_optcom_denoised.T, rcond=None)[0]
+            data_optcom_denoised[i_voxel, :] -= np.dot(
+                np.atleast_2d(betas).T,
+                design_matrix.values,
+            )
+    else:
+        betas = np.linalg.lstsq(confounds_df.values, data_optcom_denoised.T, rcond=None)[0]
+        data_optcom_denoised -= np.dot(np.atleast_2d(betas).T, confounds_df.values)
 
-    io_generator.save_file(data_oc, "has gs combined img")
-    dm_optcom = utils.unmask(tsoc_nogs, temporal_mean_mask)
-    io_generator.save_file(dm_optcom, "removed gs combined img")
+    # io_generator.save_file(data_oc, "has gs combined img")
+    data_optcom_denoised = utils.unmask(data_optcom_denoised, mask)
+    # io_generator.save_file(data_optcom_denoised, "removed regressors combined img")
 
-    # Project glbase out of each echo
-    dm_catd = data_cat.copy()  # don't overwrite data_cat
+    # Project confounds out of each echo
+    data_cat_denoised = data_cat.copy()  # don't overwrite data_cat
     for echo in range(n_echos):
-        dat = dm_catd[:, echo, :][temporal_mean_mask]
-        sol = np.linalg.lstsq(np.atleast_2d(glbase), dat.T, rcond=None)[0]
-        e_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T, np.atleast_2d(glbase.T[dtrank]))
-        dm_catd[:, echo, :] = utils.unmask(e_nogs, temporal_mean_mask)
+        echo_denoised = data_cat_denoised[:, echo, :][mask]
+        if voxelwise_confounds:
+            for i_voxel in trange(echo_denoised.shape[0], desc=f"Denoise echo {echo + 1}"):
+                design_matrix = confounds_df.copy()
+                for confound_name, confound_arr in voxelwise_confounds.items():
+                    design_matrix[confound_name] = confound_arr[i_voxel, :]
+
+                betas = np.linalg.lstsq(design_matrix.values, echo_denoised.T, rcond=None)[0]
+                echo_denoised[i_voxel, :] -= np.dot(np.atleast_2d(betas).T, design_matrix.values)
+        else:
+            betas = np.linalg.lstsq(confounds_df.values, echo_denoised.T, rcond=None)[0]
+            echo_denoised -= np.dot(np.atleast_2d(betas).T, confounds_df.values)
+
+        data_cat_denoised[:, echo, :] = utils.unmask(echo_denoised, mask)
 
-    return dm_catd, dm_optcom
+    return data_cat_denoised, data_optcom_denoised

From 39bca2037c165362f8efe776ccc6e95681a30789 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Mon, 22 Apr 2024 10:34:30 -0400
Subject: [PATCH 04/12] Handle the means.

---
 tedana/workflows/denoise_echoes.py | 45 ++++++++++++++++++++----------
 1 file changed, 31 insertions(+), 14 deletions(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index 10fbdbf93..3bc557041 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -307,6 +307,7 @@ def denoise_echoes_workflow(
         LGR.info("Computing adaptive mask")
     else:
         LGR.info("Using user-defined mask")
+
     mask, adaptive_mask = utils.make_adaptive_mask(
         data_cat,
         mask=mask,
@@ -342,21 +343,22 @@ def denoise_echoes_workflow(
     s0_full[s0_full < 0] = 0
     t2s_full[t2s_full < 0] = 0
 
-    io_generator.save_file(
-        utils.millisec2sec(t2s_full),
-        "t2star img",
-    )
+    io_generator.save_file(utils.millisec2sec(t2s_full), "t2star img")
     io_generator.save_file(s0_full, "s0 img")
-    io_generator.save_file(
-        utils.millisec2sec(t2s_limited),
-        "limited t2star img",
-    )
-    io_generator.save_file(
-        s0_limited,
-        "limited s0 img",
-    )
+    io_generator.save_file(utils.millisec2sec(t2s_limited), "limited t2star img")
+    io_generator.save_file(s0_limited, "limited s0 img")
     io_generator.save_file(data_oc, "combined img")
 
+    if confounds is not None:
+        data_cat_denoised, data_optcom_denoised = denoise_echoes(
+            data_cat=data_cat,
+            data_oc=data_oc,
+            mask=mask,
+            confounds=confounds,
+        )
+        io_generator.save_file(data_cat_denoised, "echo img")
+        io_generator.save_file(data_optcom_denoised, "optcom img")
+
     # Write out BIDS-compatible description file
     derivative_metadata = {
         "Name": "t2smap Outputs",
@@ -413,11 +415,12 @@ def _main(argv=None):
 
 
 def denoise_echoes(
+    *,
     data_cat: np.ndarray,
     data_oc: np.ndarray,
     mask: np.ndarray,
     confounds: dict[str, str],
-) -> tuple[np.ndarray, np.ndarray]:
+) -> tuple[np.ndarray, np.ndarray, pd.DataFrame]:
     """Denoise echoes using external regressors.
 
     TODO: Calculate confound-wise echo-dependence metrics.
@@ -440,6 +443,8 @@ def denoise_echoes(
         Denoised concatenated data across echoes.
     data_optcom_denoised : :obj:`numpy.ndarray` of shape (n_samples, n_volumes)
         Denoised optimally combined data.
+    metrics : :obj:`pandas.DataFrame` of shape (n_volumes, n_confounds)
+        Metrics of confound-wise echo-dependence.
     """
     LGR.info("Applying a priori confound regression")
     RepLGR.info("Confounds were regressed out of the multi-echo and optimally combined datasets.")
@@ -521,6 +526,9 @@ def denoise_echoes(
         betas = np.linalg.lstsq(confounds_df.values, data_optcom_denoised.T, rcond=None)[0]
         data_optcom_denoised -= np.dot(np.atleast_2d(betas).T, confounds_df.values)
 
+    # Add the temporal mean back
+    data_optcom_denoised += temporal_mean[mask, np.newaxis]
+
     # io_generator.save_file(data_oc, "has gs combined img")
     data_optcom_denoised = utils.unmask(data_optcom_denoised, mask)
     # io_generator.save_file(data_optcom_denoised, "removed regressors combined img")
@@ -529,6 +537,9 @@ def denoise_echoes(
     data_cat_denoised = data_cat.copy()  # don't overwrite data_cat
     for echo in range(n_echos):
         echo_denoised = data_cat_denoised[:, echo, :][mask]
+        # Remove the temporal mean
+        temporal_mean = echo_denoised.mean(axis=-1)
+        echo_denoised -= temporal_mean
         if voxelwise_confounds:
             for i_voxel in trange(echo_denoised.shape[0], desc=f"Denoise echo {echo + 1}"):
                 design_matrix = confounds_df.copy()
@@ -541,6 +552,12 @@ def denoise_echoes(
             betas = np.linalg.lstsq(confounds_df.values, echo_denoised.T, rcond=None)[0]
             echo_denoised -= np.dot(np.atleast_2d(betas).T, confounds_df.values)
 
+        # Add the temporal mean back
+        echo_denoised += temporal_mean
+
         data_cat_denoised[:, echo, :] = utils.unmask(echo_denoised, mask)
 
-    return data_cat_denoised, data_optcom_denoised
+    # TODO: Calculate metrics of confound-wise echo-dependence
+    metrics = pd.DataFrame()
+
+    return data_cat_denoised, data_optcom_denoised, metrics

From cf2be7fe1cdea0ed5bd2f188a2ab89a410a9c33a Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Tue, 30 Apr 2024 16:51:28 -0400
Subject: [PATCH 05/12] Update denoise_echoes.py

---
 tedana/workflows/denoise_echoes.py | 84 ++++++++++++++++++++++--------
 1 file changed, 63 insertions(+), 21 deletions(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index 3bc557041..d95d0de96 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -13,6 +13,7 @@
 from tqdm import trange
 
 from tedana import __version__, combine, decay, io, utils
+from tedana.metrics.collect import generate_metrics
 from tedana.workflows.parser_utils import is_valid_file
 
 LGR = logging.getLogger("GENERAL")
@@ -334,10 +335,16 @@ def denoise_echoes_workflow(
 
     LGR.info("Computing optimal combination")
     # optimally combine data
-    data_oc = combine.make_optcom(data_cat, tes, adaptive_mask, t2s=t2s_full, combmode=combmode)
+    data_optcom = combine.make_optcom(
+        data_cat,
+        tes,
+        adaptive_mask,
+        t2s=t2s_full,
+        combmode=combmode,
+    )
 
     # clean up numerical errors
-    for arr in (data_oc, s0_full, t2s_full):
+    for arr in (data_optcom, s0_full, t2s_full):
         np.nan_to_num(arr, copy=False)
 
     s0_full[s0_full < 0] = 0
@@ -347,17 +354,22 @@ def denoise_echoes_workflow(
     io_generator.save_file(s0_full, "s0 img")
     io_generator.save_file(utils.millisec2sec(t2s_limited), "limited t2star img")
     io_generator.save_file(s0_limited, "limited s0 img")
-    io_generator.save_file(data_oc, "combined img")
+    io_generator.save_file(data_optcom, "combined img")
 
     if confounds is not None:
-        data_cat_denoised, data_optcom_denoised = denoise_echoes(
+        data_cat_denoised, data_optcom_denoised, metrics = denoise_echoes(
             data_cat=data_cat,
-            data_oc=data_oc,
+            data_optcom=data_optcom,
             mask=mask,
             confounds=confounds,
         )
         io_generator.save_file(data_cat_denoised, "echo img")
         io_generator.save_file(data_optcom_denoised, "optcom img")
+        name = os.path.join(
+            io_generator.out_dir,
+            io_generator.prefix + "desc-external_metrics.tsv",
+        )
+        io_generator.save_tsv(metrics, name)
 
     # Write out BIDS-compatible description file
     derivative_metadata = {
@@ -417,9 +429,11 @@ def _main(argv=None):
 def denoise_echoes(
     *,
     data_cat: np.ndarray,
-    data_oc: np.ndarray,
-    mask: np.ndarray,
+    tes: np.ndarray,
+    data_optcom: np.ndarray,
+    adaptive_mask: np.ndarray,
     confounds: dict[str, str],
+    io_generator: io.OutputGenerator,
 ) -> tuple[np.ndarray, np.ndarray, pd.DataFrame]:
     """Denoise echoes using external regressors.
 
@@ -429,13 +443,17 @@ def denoise_echoes(
     ----------
     data_cat : :obj:`numpy.ndarray` of shape (n_samples, n_echos, n_volumes)
         Concatenated data across echoes.
-    data_oc : :obj:`numpy.ndarray` of shape (n_samples, n_volumes)
+    tes : :obj:`numpy.ndarray` of shape (n_echos,)
+        Echo times in milliseconds.
+    data_optcom : :obj:`numpy.ndarray` of shape (n_samples, n_volumes)
         Optimally combined data across echoes.
-    mask : :obj:`numpy.ndarray` of shape (n_samples,)
-        Binary mask of voxels to include in TE Dependent ANAlysis.
+    adaptive_mask : :obj:`numpy.ndarray` of shape (n_samples,)
+        Adaptive mask of voxels to include in TE Dependent ANAlysis.
     confounds : :obj:`dict`
         Files defining confounds to regress from the echo-wise data.
         Keys indicate the names to use for the confounds.
+    io_generator : :obj:`~tedana.io.OutputGenerator`
+        Output generator.
 
     Returns
     -------
@@ -449,15 +467,16 @@ def denoise_echoes(
     LGR.info("Applying a priori confound regression")
     RepLGR.info("Confounds were regressed out of the multi-echo and optimally combined datasets.")
     n_samples, n_echos, n_volumes = data_cat.shape
-    if data_cat.shape[0] != data_oc.shape[0]:
+    mask = adaptive_mask >= 1
+    if data_cat.shape[0] != data_optcom.shape[0]:
         raise ValueError(
-            f"First dimensions of data_cat ({data_cat.shape[0]}) and data_oc ({data_oc.shape[0]}) "
+            f"First dimensions of data_cat ({data_cat.shape[0]}) and data_optcom ({data_optcom.shape[0]}) "
             "do not match"
         )
-    elif data_cat.shape[2] != data_oc.shape[1]:
+    elif data_cat.shape[2] != data_optcom.shape[1]:
         raise ValueError(
             f"Third dimension of data_cat ({data_cat.shape[2]}) does not match second dimension "
-            f"of data_oc ({data_oc.shape[1]})"
+            f"of data_optcom ({data_optcom.shape[1]})"
         )
 
     confound_dfs = []
@@ -481,7 +500,7 @@ def denoise_echoes(
                 # Find the time course for the spatial map
                 confound_timeseries = np.linalg.lstsq(
                     np.atleast_2d(confound_arr).T,
-                    data_oc,
+                    data_optcom,
                     rcond=None,
                 )[0]
                 confound_timeseries = stats.zscore(confound_timeseries, axis=None)
@@ -505,12 +524,38 @@ def denoise_echoes(
     if confound_dfs:
         confounds_df = pd.concat(confound_dfs, axis=1)
         confounds_df = confounds_df.fillna(0)
+
+        # Calculate dependence metrics for non-voxel-wise confounds
+        # TODO: Support metrics for voxel-wise confounds too
+        required_metrics = [
+            "kappa",
+            "rho",
+            "countnoise",
+            "countsigFT2",
+            "countsigFS0",
+            "dice_FT2",
+            "dice_FS0",
+            "signal-noise_t",
+            "variance explained",
+            "normalized variance explained",
+            "d_table_score",
+        ]
+        metrics = generate_metrics(
+            data_cat=data_cat,
+            data_optcom=data_optcom,
+            mixing=confounds_df.values,
+            adaptive_mask=adaptive_mask,
+            tes=tes,
+            io_generator=io_generator,
+            label="external",
+            metrics=required_metrics,
+        )
     else:
         confounds_df = pd.DataFrame(index=range(n_volumes))
 
     # Project confounds out of optimally combined data
-    temporal_mean = data_oc.mean(axis=-1)  # temporal mean
-    data_optcom_denoised = data_oc[mask] - temporal_mean[mask, np.newaxis]
+    temporal_mean = data_optcom.mean(axis=-1)  # temporal mean
+    data_optcom_denoised = data_optcom[mask] - temporal_mean[mask, np.newaxis]
     if voxelwise_confounds:
         for i_voxel in trange(data_optcom_denoised.shape[0], desc="Denoise optimally combined"):
             design_matrix = confounds_df.copy()
@@ -529,7 +574,7 @@ def denoise_echoes(
     # Add the temporal mean back
     data_optcom_denoised += temporal_mean[mask, np.newaxis]
 
-    # io_generator.save_file(data_oc, "has gs combined img")
+    # io_generator.save_file(data_optcom, "has gs combined img")
     data_optcom_denoised = utils.unmask(data_optcom_denoised, mask)
     # io_generator.save_file(data_optcom_denoised, "removed regressors combined img")
 
@@ -557,7 +602,4 @@ def denoise_echoes(
 
         data_cat_denoised[:, echo, :] = utils.unmask(echo_denoised, mask)
 
-    # TODO: Calculate metrics of confound-wise echo-dependence
-    metrics = pd.DataFrame()
-
     return data_cat_denoised, data_optcom_denoised, metrics

From a5fa0e7ef5047a721ce00342e08427d3e4cd75f5 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Sun, 5 May 2024 17:02:48 -0400
Subject: [PATCH 06/12] More work.

---
 pyproject.toml                     |  1 +
 tedana/workflows/__init__.py       |  8 ++++++-
 tedana/workflows/denoise_echoes.py | 36 +++++++++++++++++++++++-------
 3 files changed, 36 insertions(+), 9 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index f757350e7..806b527ff 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -76,6 +76,7 @@ all = ["tedana[dev,doc,tests]"]
 ica_reclassify = "tedana.workflows.ica_reclassify:_main"
 t2smap = "tedana.workflows.t2smap:_main"
 tedana = "tedana.workflows.tedana:_main"
+denoise_echoes = "tedana.workflows.denoise_echoes:_main"
 
 #
 # Hatch configurations
diff --git a/tedana/workflows/__init__.py b/tedana/workflows/__init__.py
index 04b284826..6dbb5a7fd 100644
--- a/tedana/workflows/__init__.py
+++ b/tedana/workflows/__init__.py
@@ -1,8 +1,14 @@
 # emacs: -*- mode: python-mode; py-indent-offset: 4; tab-width: 4; indent-tabs-mode: nil -*-
 # ex: set sts=4 ts=4 sw=4 et:
 """Command line interfaces and workflows."""
+from tedana.workflows.denoise_echoes import denoise_echoes_workflow
 from tedana.workflows.ica_reclassify import ica_reclassify_workflow
 from tedana.workflows.t2smap import t2smap_workflow
 from tedana.workflows.tedana import tedana_workflow
 
-__all__ = ["tedana_workflow", "t2smap_workflow", "ica_reclassify_workflow"]
+__all__ = [
+    "tedana_workflow",
+    "t2smap_workflow",
+    "ica_reclassify_workflow",
+    "denoise_echoes_workflow",
+]
diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index d95d0de96..641eeeb01 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -5,6 +5,7 @@
 import os
 import os.path as op
 import sys
+from pathlib import Path
 
 import numpy as np
 import pandas as pd
@@ -28,6 +29,24 @@ def _get_parser():
     parser.parse_args() : argparse dict
     """
     parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+
+    class ToDict(argparse.Action):
+        def __call__(self, parser, namespace, values, option_string=None):
+            d = {}
+            for spec in values:
+                try:
+                    name, loc = spec.split("=")
+                    loc = Path(loc)
+                except ValueError:
+                    loc = Path(spec)
+                    name = loc.name
+
+                if name in d:
+                    raise ValueError(f"Received duplicate derivative name: {name}")
+
+                d[name] = loc
+            setattr(namespace, self.dest, d)
+
     # Argument parser follow template provided by RalphyZ
     # https://stackoverflow.com/a/43456577
     optional = parser._action_groups.pop()
@@ -59,9 +78,9 @@ def _get_parser():
     required.add_argument(
         "--confounds",
         dest="confounds",
-        nargs="+",
-        metavar="FILE",
+        action=ToDict,
         type=lambda x: is_valid_file(parser, x),
+        nargs="+",
         help="Files defining confounds to regress from the echo-wise data.",
     )
     optional.add_argument(
@@ -182,8 +201,7 @@ def denoise_echoes_workflow(
     quiet=False,
     t2smap_command=None,
 ):
-    """
-    Estimate T2 and S0, and optimally combine data across TEs.
+    """Estimate T2 and S0, and optimally combine data across TEs.
 
     Please remember to cite :footcite:t:`dupre2021te`.
 
@@ -196,7 +214,7 @@ def denoise_echoes_workflow(
         List of echo times associated with data in milliseconds.
     out_dir : :obj:`str`, optional
         Output directory.
-    confounds : :obj:`list` of :obj:`str`
+    confounds : :obj:`dict` of :obj:`str`
         Files defining confounds to regress from the echo-wise data.
     mask : :obj:`str`, optional
         Binary mask of voxels to include in TE Dependent ANAlysis. Must be spatially
@@ -359,9 +377,11 @@ def denoise_echoes_workflow(
     if confounds is not None:
         data_cat_denoised, data_optcom_denoised, metrics = denoise_echoes(
             data_cat=data_cat,
+            tes=tes,
             data_optcom=data_optcom,
-            mask=mask,
+            adaptive_mask=adaptive_mask,
             confounds=confounds,
+            io_generator=io_generator,
         )
         io_generator.save_file(data_cat_denoised, "echo img")
         io_generator.save_file(data_optcom_denoised, "optcom img")
@@ -470,8 +490,8 @@ def denoise_echoes(
     mask = adaptive_mask >= 1
     if data_cat.shape[0] != data_optcom.shape[0]:
         raise ValueError(
-            f"First dimensions of data_cat ({data_cat.shape[0]}) and data_optcom ({data_optcom.shape[0]}) "
-            "do not match"
+            f"First dimensions of data_cat ({data_cat.shape[0]}) and "
+            f"data_optcom ({data_optcom.shape[0]}) do not match"
         )
     elif data_cat.shape[2] != data_optcom.shape[1]:
         raise ValueError(

From 2cbf0727b627fb664c0ffd9e7c772726d6f37ed1 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Sun, 5 May 2024 18:54:03 -0400
Subject: [PATCH 07/12] More work.

---
 tedana/metrics/collect.py            | 14 ++++++++++----
 tedana/resources/config/outputs.json |  6 +++---
 tedana/workflows/denoise_echoes.py   | 26 ++++++++++++++++----------
 3 files changed, 29 insertions(+), 17 deletions(-)

diff --git a/tedana/metrics/collect.py b/tedana/metrics/collect.py
index e71e5ee4a..cff5c7ace 100644
--- a/tedana/metrics/collect.py
+++ b/tedana/metrics/collect.py
@@ -82,6 +82,9 @@ def generate_metrics(
             "match."
         )
 
+    if isinstance(label, list):
+        assert len(label) == mixing.shape[1]
+
     # Derive mask from thresholded adaptive mask
     mask = adaptive_mask >= 3
 
@@ -109,10 +112,13 @@ def generate_metrics(
     # throughout this function
     n_components = mixing.shape[1]
     comptable = pd.DataFrame(index=np.arange(n_components, dtype=int))
-    comptable["Component"] = [
-        io.add_decomp_prefix(comp, prefix=label, max_value=comptable.shape[0])
-        for comp in comptable.index.values
-    ]
+    if isinstance(label, str):
+        comptable["Component"] = [
+            io.add_decomp_prefix(comp, prefix=label, max_value=comptable.shape[0])
+            for comp in comptable.index.values
+        ]
+    else:
+        comptable["Component"] = label
 
     # Metric maps
     # Maps will be stored as arrays in an easily-indexable dictionary
diff --git a/tedana/resources/config/outputs.json b/tedana/resources/config/outputs.json
index 0ada7b33d..fccda12f0 100644
--- a/tedana/resources/config/outputs.json
+++ b/tedana/resources/config/outputs.json
@@ -121,15 +121,15 @@
     },
     "high kappa ts split img": {
         "orig": "hik_ts_e{echo}",
-        "bidsv1.5.0": "echo-{echo}_desc-Accepted_bold"
+        "bidsv1.5.0": "echo-{echo}_desc-accepted_bold"
     },
     "low kappa ts split img": {
         "orig": "lowk_ts_e{echo}",
-        "bidsv1.5.0": "echo-{echo}_desc-Rejected_bold"
+        "bidsv1.5.0": "echo-{echo}_desc-rejected_bold"
     },
     "denoised ts split img": {
         "orig": "dn_ts_e{echo}",
-        "bidsv1.5.0": "echo-{echo}_desc-Denoised_bold"
+        "bidsv1.5.0": "echo-{echo}_desc-denoised_bold"
     },
     "gs img": {
         "orig": "T1gs",
diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index 641eeeb01..0ad8e7961 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -383,8 +383,14 @@ def denoise_echoes_workflow(
             confounds=confounds,
             io_generator=io_generator,
         )
-        io_generator.save_file(data_cat_denoised, "echo img")
-        io_generator.save_file(data_optcom_denoised, "optcom img")
+        io_generator.save_file(data_optcom_denoised, "denoised ts img")
+        for echo in range(n_echos):
+            io_generator.save_file(
+                data_cat_denoised[:, echo, :],
+                "denoised ts split img",
+                echo=echo + 1,
+            )
+
         name = os.path.join(
             io_generator.out_dir,
             io_generator.prefix + "desc-external_metrics.tsv",
@@ -393,7 +399,7 @@ def denoise_echoes_workflow(
 
     # Write out BIDS-compatible description file
     derivative_metadata = {
-        "Name": "t2smap Outputs",
+        "Name": "denoise_echoes Outputs",
         "BIDSVersion": "1.5.0",
         "DatasetType": "derivative",
         "GeneratedBy": [
@@ -567,15 +573,15 @@ def denoise_echoes(
             adaptive_mask=adaptive_mask,
             tes=tes,
             io_generator=io_generator,
-            label="external",
+            label=confounds_df.columns.tolist(),
             metrics=required_metrics,
         )
     else:
         confounds_df = pd.DataFrame(index=range(n_volumes))
 
     # Project confounds out of optimally combined data
-    temporal_mean = data_optcom.mean(axis=-1)  # temporal mean
-    data_optcom_denoised = data_optcom[mask] - temporal_mean[mask, np.newaxis]
+    temporal_mean = data_optcom.mean(axis=-1, keepdims=True)  # temporal mean
+    data_optcom_denoised = data_optcom[mask] - temporal_mean[mask]
     if voxelwise_confounds:
         for i_voxel in trange(data_optcom_denoised.shape[0], desc="Denoise optimally combined"):
             design_matrix = confounds_df.copy()
@@ -589,10 +595,10 @@ def denoise_echoes(
             )
     else:
         betas = np.linalg.lstsq(confounds_df.values, data_optcom_denoised.T, rcond=None)[0]
-        data_optcom_denoised -= np.dot(np.atleast_2d(betas).T, confounds_df.values)
+        data_optcom_denoised -= np.dot(np.atleast_2d(betas).T, confounds_df.values.T)
 
     # Add the temporal mean back
-    data_optcom_denoised += temporal_mean[mask, np.newaxis]
+    data_optcom_denoised += temporal_mean[mask]
 
     # io_generator.save_file(data_optcom, "has gs combined img")
     data_optcom_denoised = utils.unmask(data_optcom_denoised, mask)
@@ -603,7 +609,7 @@ def denoise_echoes(
     for echo in range(n_echos):
         echo_denoised = data_cat_denoised[:, echo, :][mask]
         # Remove the temporal mean
-        temporal_mean = echo_denoised.mean(axis=-1)
+        temporal_mean = echo_denoised.mean(axis=-1, keepdims=True)
         echo_denoised -= temporal_mean
         if voxelwise_confounds:
             for i_voxel in trange(echo_denoised.shape[0], desc=f"Denoise echo {echo + 1}"):
@@ -615,7 +621,7 @@ def denoise_echoes(
                 echo_denoised[i_voxel, :] -= np.dot(np.atleast_2d(betas).T, design_matrix.values)
         else:
             betas = np.linalg.lstsq(confounds_df.values, echo_denoised.T, rcond=None)[0]
-            echo_denoised -= np.dot(np.atleast_2d(betas).T, confounds_df.values)
+            echo_denoised -= np.dot(np.atleast_2d(betas).T, confounds_df.values.T)
 
         # Add the temporal mean back
         echo_denoised += temporal_mean

From 2ec5c74dd47a08607d521803844874c4268bea5a Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Sun, 19 May 2024 10:15:25 -0400
Subject: [PATCH 08/12] Update docstrings.

---
 tedana/workflows/denoise_echoes.py | 35 +++++++++++++++---------------
 1 file changed, 17 insertions(+), 18 deletions(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index 0ad8e7961..f84d4ac29 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -199,9 +199,9 @@ def denoise_echoes_workflow(
     combmode="t2s",
     debug=False,
     quiet=False,
-    t2smap_command=None,
+    denoise_echoes_command=None,
 ):
-    """Estimate T2 and S0, and optimally combine data across TEs.
+    """Denoise individual echoes using external regressors.
 
     Please remember to cite :footcite:t:`dupre2021te`.
 
@@ -219,6 +219,8 @@ def denoise_echoes_workflow(
     mask : :obj:`str`, optional
         Binary mask of voxels to include in TE Dependent ANAlysis. Must be spatially
         aligned with `data`.
+    prefix
+    convention
     masktype : :obj:`list` with 'dropout' and/or 'decay' or None, optional
         Method(s) by which to define the adaptive mask. Default is ["dropout"].
     fittype : {'loglin', 'curvefit'}, optional
@@ -234,7 +236,7 @@ def denoise_echoes_workflow(
         Default is 'all'.
     combmode : {'t2s', 'paid'}, optional
         Combination scheme for TEs: 't2s' (Posse 1999, default), 'paid' (Poser).
-    t2smap_command : :obj:`str`, optional
+    denoise_echoes_command : :obj:`str`, optional
         The command used to run t2smap. Default is None.
 
     Other Parameters
@@ -286,20 +288,20 @@ def denoise_echoes_workflow(
     LGR.info(f"Using output directory: {out_dir}")
 
     # Save command into sh file, if the command-line interface was used
-    if t2smap_command is not None:
-        command_file = open(os.path.join(out_dir, "t2smap_call.sh"), "w")
-        command_file.write(t2smap_command)
+    if denoise_echoes_command is not None:
+        command_file = open(os.path.join(out_dir, "denoise_echoes_call.sh"), "w")
+        command_file.write(denoise_echoes_command)
         command_file.close()
     else:
         # Get variables passed to function if the tedana command is None
         variables = ", ".join(f"{name}={value}" for name, value in locals().items())
         # From variables, remove everything after ", tedana_command"
-        variables = variables.split(", t2smap_command")[0]
-        t2smap_command = f"t2smap_workflow({variables})"
+        variables = variables.split(", denoise_echoes_command")[0]
+        denoise_echoes_command = f"denoise_echoes_workflow({variables})"
 
     # Save system info to json
     info_dict = utils.get_system_version_info()
-    info_dict["Command"] = t2smap_command
+    info_dict["Command"] = denoise_echoes_command
 
     # ensure tes are in appropriate format
     tes = [float(te) for te in tes]
@@ -406,10 +408,7 @@ def denoise_echoes_workflow(
             {
                 "Name": "t2smap",
                 "Version": __version__,
-                "Description": (
-                    "A pipeline estimating T2* from multi-echo fMRI data and "
-                    "combining data across echoes."
-                ),
+                "Description": "A pipeline denoising echo-wise data using external regressors.",
                 "CodeURL": "https://github.com/ME-ICA/tedana",
                 "Node": {
                     "Name": info_dict["Node"],
@@ -434,18 +433,18 @@ def denoise_echoes_workflow(
 
 
 def _main(argv=None):
-    """Run the t2smap workflow."""
+    """Run the denoise_echoes workflow."""
     if argv:
         # relevant for tests when CLI called with t2smap_cli._main(args)
-        t2smap_command = "t2smap " + " ".join(argv)
+        denoise_echoes_command = "denoise_echoes " + " ".join(argv)
     else:
-        t2smap_command = "t2smap " + " ".join(sys.argv[1:])
+        denoise_echoes_command = "denoise_echoes " + " ".join(sys.argv[1:])
     options = _get_parser().parse_args(argv)
     kwargs = vars(options)
     n_threads = kwargs.pop("n_threads")
     n_threads = None if n_threads == -1 else n_threads
     with threadpool_limits(limits=n_threads, user_api=None):
-        denoise_echoes_workflow(**kwargs, t2smap_command=t2smap_command)
+        denoise_echoes_workflow(**kwargs, denoise_echoes_command=denoise_echoes_command)
 
 
 if __name__ == "__main__":
@@ -463,7 +462,7 @@ def denoise_echoes(
 ) -> tuple[np.ndarray, np.ndarray, pd.DataFrame]:
     """Denoise echoes using external regressors.
 
-    TODO: Calculate confound-wise echo-dependence metrics.
+    TODO: Support voxel-wise external regressors.
 
     Parameters
     ----------

From a7230ac39628c21d8ec16f98ef0b7909a5a914d2 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Sun, 19 May 2024 10:38:30 -0400
Subject: [PATCH 09/12] Update denoise_echoes.py

---
 tedana/workflows/denoise_echoes.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index f84d4ac29..726b8ac08 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -6,6 +6,7 @@
 import os.path as op
 import sys
 from pathlib import Path
+from typing import Dict
 
 import numpy as np
 import pandas as pd
@@ -457,7 +458,7 @@ def denoise_echoes(
     tes: np.ndarray,
     data_optcom: np.ndarray,
     adaptive_mask: np.ndarray,
-    confounds: dict[str, str],
+    confounds: Dict[str, str],
     io_generator: io.OutputGenerator,
 ) -> tuple[np.ndarray, np.ndarray, pd.DataFrame]:
     """Denoise echoes using external regressors.

From abc5168a71cf42ac19c9d632af12d616d046c0ac Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Sun, 19 May 2024 10:39:26 -0400
Subject: [PATCH 10/12] Update denoise_echoes.py

---
 tedana/workflows/denoise_echoes.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index 726b8ac08..70e5c48e4 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -436,7 +436,7 @@ def denoise_echoes_workflow(
 def _main(argv=None):
     """Run the denoise_echoes workflow."""
     if argv:
-        # relevant for tests when CLI called with t2smap_cli._main(args)
+        # relevant for tests when CLI called with denoise_echoes_cli._main(args)
         denoise_echoes_command = "denoise_echoes " + " ".join(argv)
     else:
         denoise_echoes_command = "denoise_echoes " + " ".join(sys.argv[1:])

From 0afb225c7378a8febb699f4684d0840f035220c3 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Sun, 19 May 2024 10:40:08 -0400
Subject: [PATCH 11/12] Update denoise_echoes.py

---
 tedana/workflows/denoise_echoes.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index 70e5c48e4..2edc08ab1 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -6,7 +6,7 @@
 import os.path as op
 import sys
 from pathlib import Path
-from typing import Dict
+from typing import Dict, Tuple
 
 import numpy as np
 import pandas as pd
@@ -460,7 +460,7 @@ def denoise_echoes(
     adaptive_mask: np.ndarray,
     confounds: Dict[str, str],
     io_generator: io.OutputGenerator,
-) -> tuple[np.ndarray, np.ndarray, pd.DataFrame]:
+) -> Tuple[np.ndarray, np.ndarray, pd.DataFrame]:
     """Denoise echoes using external regressors.
 
     TODO: Support voxel-wise external regressors.

From f9785ae76ecf96aa21f3ef4d984a86a38483fd08 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Sun, 19 May 2024 10:44:41 -0400
Subject: [PATCH 12/12] Update denoise_echoes.py

---
 tedana/workflows/denoise_echoes.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tedana/workflows/denoise_echoes.py b/tedana/workflows/denoise_echoes.py
index 2edc08ab1..d140bfb82 100644
--- a/tedana/workflows/denoise_echoes.py
+++ b/tedana/workflows/denoise_echoes.py
@@ -32,7 +32,7 @@ def _get_parser():
     parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
     class ToDict(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
+        def __call__(self, parser, namespace, values, option_string=None):  # noqa: U100
             d = {}
             for spec in values:
                 try: