[ENH] Add ability to re-run ICA when no BOLD components are found

tedana/decomposition/ica.py

@@ -38,6 +38,8 @@ def tedica(data, n_components, fixed_seed, maxit=500, maxrestart=10):
     mmix : (T x C) :obj:`numpy.ndarray`
         Z-scored mixing matrix for converting input data to component space,
         where `C` is components and `T` is the same as in `data`
+    fixed_seed : :obj:`int`
+        Random seed from final decomposition.
 
     Notes
     -----
@@ -64,16 +66,16 @@ def tedica(data, n_components, fixed_seed, maxit=500, maxrestart=10):
 
             w = list(filter(lambda i: issubclass(i.category, UserWarning), w))
             if len(w):
-                LGR.warning('ICA attempt {0} failed to converge after {1} '
-                            'iterations'.format(i_attempt + 1, ica.n_iter_))
+                LGR.warning('ICA with random seed {0} failed to converge after {1} '
+                            'iterations'.format(fixed_seed, ica.n_iter_))
                 if i_attempt < maxrestart - 1:
                     fixed_seed += 1
                     LGR.warning('Random seed updated to {0}'.format(fixed_seed))
             else:
-                LGR.info('ICA attempt {0} converged in {1} '
-                         'iterations'.format(i_attempt + 1, ica.n_iter_))
+                LGR.info('ICA with random seed {0} converged in {1} '
+                         'iterations'.format(fixed_seed, ica.n_iter_))
                 break
 
     mmix = ica.mixing_
     mmix = stats.zscore(mmix, axis=0)
-    return mmix
+    return mmix, fixed_seed

tedana/tests/data/fiu_four_echo_outputs.txt

@@ -1,15 +1,22 @@
 T1gs.nii.gz
 adaptive_mask.nii.gz
 betas_OC.nii.gz
+betas_hik_OC.nii.gz
 betas_hik_OC_MIR.nii.gz
 dn_ts_OC.nii.gz
 dn_ts_OC_MIR.nii.gz
 dn_ts_e1.nii.gz
 dn_ts_e2.nii.gz
 dn_ts_e3.nii.gz
 dn_ts_e4.nii.gz
+feats_OC2.nii.gz
 glsig.1D
+hik_ts_OC.nii.gz
 hik_ts_OC_MIR.nii.gz
+hik_ts_e1.nii.gz
+hik_ts_e2.nii.gz
+hik_ts_e3.nii.gz
+hik_ts_e4.nii.gz
 ica_components.nii.gz
 ica_decomposition.json
 ica_mixing.tsv

tedana/workflows/tedana.py

@@ -514,21 +514,45 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
                                                 out_dir=out_dir,
                                                 verbose=verbose,
                                                 low_mem=low_mem)
-        mmix_orig = decomposition.tedica(dd, n_components, fixed_seed,
-                                         maxit, maxrestart)
-
         if verbose:
             io.filewrite(utils.unmask(dd, mask),
                          op.join(out_dir, 'ts_OC_whitened.nii.gz'), 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)
-        comptable, metric_maps, betas, mmix = metrics.dependence_metrics(
-                    catd, data_oc, mmix_orig, masksum, tes,
-                    ref_img, reindex=True, label='meica_', out_dir=out_dir,
-                    algorithm='kundu_v2', verbose=verbose)
+        # Perform ICA, calculate metrics, and apply decision tree
+        # Restart when ICA fails to converge or too few BOLD components found
+        keep_restarting = True
+        n_restarts = 0
+        seed = fixed_seed
+        while keep_restarting:
+            mmix_orig, seed = decomposition.tedica(
+                dd, n_components, seed,
+                maxit, maxrestart=(maxrestart - n_restarts)
+            )
+            seed += 1
+            n_restarts = seed - fixed_seed
+
+            # Estimate betas and compute selection metrics for mixing matrix
+            # generated from dimensionally reduced data using full data (i.e., data
+            # with thermal noise)
+            LGR.info('Making second component selection guess from ICA results')
+            comptable, metric_maps, betas, mmix = metrics.dependence_metrics(
+                catd, data_oc, mmix_orig, masksum, tes,
+                ref_img, reindex=True, label='meica_', out_dir=out_dir,
+                algorithm='kundu_v2', verbose=verbose
+            )
+            comptable = metrics.kundu_metrics(comptable, metric_maps)
+            comptable = selection.kundu_selection_v2(comptable, n_echos, n_vols)
+
+            n_bold_comps = comptable[comptable.classification == 'accepted'].shape[0]
+            if (n_restarts < maxrestart) and (n_bold_comps == 0):
+                LGR.warning("No BOLD components found. Re-attempting ICA.")
+            elif (n_bold_comps == 0):
+                LGR.warning("No BOLD components found, but maximum number of restarts reached.")
+                keep_restarting = False
+            else:
+                keep_restarting = False
+
+        # Write out ICA files.
         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)
@@ -537,9 +561,6 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
         io.filewrite(betas_oc,
                      op.join(out_dir, 'ica_components.nii.gz'),
                      ref_img)
-
-        comptable = metrics.kundu_metrics(comptable, metric_maps)
-        comptable = selection.kundu_selection_v2(comptable, n_echos, n_vols)
     else:
         LGR.info('Using supplied mixing matrix from ICA')
         mmix_orig = pd.read_table(op.join(out_dir, 'ica_mixing.tsv')).values
@@ -561,7 +582,7 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
                      op.join(out_dir, 'ica_components.nii.gz'),
                      ref_img)
 
-    # Save decomposition
+    # Save component table
     comptable['Description'] = 'ICA fit to dimensionally-reduced optimally combined data.'
     mmix_dict = {}
     mmix_dict['Method'] = ('Independent components analysis with FastICA '