[ENH] Add verbose outputs for pipeline walkthrough

.circleci/tedana_outputs_verbose.txt

@@ -23,8 +23,16 @@ lowk_ts_e2.nii
 lowk_ts_e3.nii
 lowk_ts_e4.nii
 lowk_ts_e5.nii
+meica_R2_pred.nii
+meica_S0_pred.nii
+meica_betas_catd.nii
+meica_metric_weights.nii
 meica_mix.1D
 mepca_OC_components.nii
+mepca_R2_pred.nii
+mepca_S0_pred.nii
+mepca_betas_catd.nii
+mepca_metric_weights.nii
 mepca_mix.1D
 pcastate.pkl
 s0v.nii
@@ -34,3 +42,4 @@ t2ss.nii
 t2sv.nii
 t2svG.nii
 ts_OC.nii
+tsoc_whitened.nii

tedana/decomposition/eigendecomp.py

@@ -165,7 +165,8 @@ def kundu_tedpca(ct_df, n_echos, kdaw, rdaw, stabilize=False):
 
 
 def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG,
-           ref_img, tes, method='mle', ste=0, kdaw=10., rdaw=1., wvpca=False):
+           ref_img, tes, method='mle', ste=0, kdaw=10., rdaw=1., wvpca=False,
+           verbose=False):
     """
     Use principal components analysis (PCA) to identify and remove thermal
     noise from multi-echo data.
@@ -199,6 +200,8 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG,
         a subset of echos. Default: 0
     wvpca : :obj:`bool`, optional
         Whether to apply wavelet denoising to data. Default: False
+    verbose : :obj:`bool`, optional
+        Whether to output files from fitmodels_direct or not. Default: False
 
     Returns
     -------
@@ -322,7 +325,8 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG,
         LGR.info('Making initial component selection guess from PCA results')
         _, ct_df, betasv, v_T = model.fitmodels_direct(
                     catd, comp_ts.T, eimum, t2s, t2sG, tes, combmode, ref_img,
-                    mmixN=vTmixN, full_sel=False)
+                    mmixN=vTmixN, full_sel=False, label='mepca_',
+                    verbose=verbose)
         # varex_norm overrides normalized varex computed by fitmodels_direct
         ct_df['normalized variance explained'] = varex_norm
 

tedana/model/fit.py

@@ -2,6 +2,7 @@
 Fit models.
 """
 import logging
+import os.path as op
 
 import numpy as np
 import pandas as pd
@@ -20,7 +21,8 @@
 
 
 def fitmodels_direct(catd, mmix, mask, t2s, t2s_full, tes, combmode, ref_img,
-                     reindex=False, mmixN=None, full_sel=True):
+                     reindex=False, mmixN=None, full_sel=True, label=None,
+                     out_dir='.', verbose=False):
     """
     Fit TE-dependence and -independence models to components.
 
@@ -150,6 +152,8 @@ def fitmodels_direct(catd, mmix, mask, t2s, t2s_full, tes, combmode, ref_img,
     F_S0_clmaps = np.zeros([n_data_voxels, n_components])
     Br_R2_clmaps = np.zeros([n_voxels, n_components])
     Br_S0_clmaps = np.zeros([n_voxels, n_components])
+    pred_R2_maps = np.zeros([n_data_voxels, n_echos, n_components])
+    pred_S0_maps = np.zeros([n_data_voxels, n_echos, n_components])
 
     LGR.info('Fitting TE- and S0-dependent models to components')
     for i_comp in range(n_components):
@@ -164,6 +168,7 @@ def fitmodels_direct(catd, mmix, mask, t2s, t2s_full, tes, combmode, ref_img,
         # (S,) model coefficient map
         coeffs_S0 = (B * X1).sum(axis=0) / (X1**2).sum(axis=0)
         pred_S0 = X1 * np.tile(coeffs_S0, (n_echos, 1))
+        pred_S0_maps[:, :, i_comp] = pred_S0.T
         SSE_S0 = (B - pred_S0)**2
         SSE_S0 = SSE_S0.sum(axis=0)  # (S,) prediction error map
         F_S0 = (alpha - SSE_S0) * (n_echos - 1) / (SSE_S0)
@@ -172,6 +177,7 @@ def fitmodels_direct(catd, mmix, mask, t2s, t2s_full, tes, combmode, ref_img,
         # R2 Model
         coeffs_R2 = (B * X2).sum(axis=0) / (X2**2).sum(axis=0)
         pred_R2 = X2 * np.tile(coeffs_R2, (n_echos, 1))
+        pred_R2_maps[:, :, i_comp] = pred_R2.T
         SSE_R2 = (B - pred_R2)**2
         SSE_R2 = SSE_R2.sum(axis=0)
         F_R2 = (alpha - SSE_R2) * (n_echos - 1) / (SSE_R2)
@@ -198,6 +204,9 @@ def fitmodels_direct(catd, mmix, mask, t2s, t2s_full, tes, combmode, ref_img,
         sort_idx = comptab[:, 0].argsort()[::-1]
         comptab = comptab[sort_idx, :]
         mmix_new = mmix[:, sort_idx]
+        betas = betas[..., sort_idx]
+        pred_R2_maps = pred_R2_maps[:, :, sort_idx]
+        pred_S0_maps = pred_S0_maps[:, :, sort_idx]
         F_S0_maps = F_S0_maps[:, sort_idx]
         F_R2_maps = F_R2_maps[:, sort_idx]
         Z_maps = Z_maps[:, sort_idx]
@@ -208,6 +217,19 @@ def fitmodels_direct(catd, mmix, mask, t2s, t2s_full, tes, combmode, ref_img,
     else:
         mmix_new = mmix
 
+    if verbose:
+        # Echo-specific weight maps for each of the ICA components.
+        io.filewrite(betas, op.join(out_dir, label+'betas_catd.nii'), ref_img)
+        # Echo-specific maps of predicted values for R2 and S0 models for each
+        # component.
+        io.filewrite(utils.unmask(pred_R2_maps, t2s != 0),
+                     op.join(out_dir, label+'R2_pred.nii'), ref_img)
+        io.filewrite(utils.unmask(pred_S0_maps, t2s != 0),
+                     op.join(out_dir, label+'S0_pred.nii'), ref_img)
+        # Weight maps used to average metrics across voxels
+        io.filewrite(utils.unmask(Z_maps ** 2., mask),
+                     op.join(out_dir, label+'metric_weights.nii'), ref_img)
+
     comptab = pd.DataFrame(comptab,
                            columns=['kappa', 'rho',
                                     'variance explained',

tedana/workflows/tedana.py

@@ -302,20 +302,25 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None,
         n_components, dd = decomposition.tedpca(catd, data_oc, combmode, mask,
                                                 t2s, t2sG, ref_img,
                                                 tes=tes, method=tedpca, ste=ste,
-                                                kdaw=10., rdaw=1., wvpca=wvpca)
+                                                kdaw=10., rdaw=1., wvpca=wvpca,
+                                                verbose=verbose)
         mmix_orig, fixed_seed = decomposition.tedica(n_components, dd,
                                                      fixed_seed)
 
         if verbose:
             np.savetxt(op.join(out_dir, '__meica_mix.1D'), mmix_orig)
+            if ste == -1:
+                io.filewrite(utils.unmask(dd, mask),
+                             op.join(out_dir, 'tsoc_whitened.nii'), ref_img)
 
         LGR.info('Making second component selection guess from ICA results')
         # Estimate betas and compute selection metrics for mixing matrix
         # generated from dimensionally reduced data using full data (i.e., data
         # with thermal noise)
         seldict, comptable, betas, mmix = model.fitmodels_direct(
                     catd, mmix_orig, mask, t2s, t2sG, tes, combmode,
-                    ref_img, reindex=True)
+                    ref_img, reindex=True, label='meica_', out_dir=out_dir,
+                    verbose=verbose)
         np.savetxt(op.join(out_dir, 'meica_mix.1D'), mmix)
 
         comptable = selection.selcomps(seldict, comptable, mmix, manacc,
@@ -325,7 +330,8 @@ def tedana_workflow(data, tes, mask=None, mixm=None, ctab=None, manacc=None,
         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)
+                    ref_img, label='meica_', out_dir=out_dir,
+                    verbose=verbose)
         if ctab is None:
             comptable = selection.selcomps(seldict, comptable, mmix, manacc,
                                            n_echos)