Merge pull request #1443 from effigies/fix/compcor_nans

[MRG] FIX: Output NaN columns for CompCor on failure

Author: effigies

docs/environment.yml

@@ -20,7 +20,7 @@ dependencies:
 - python-dateutil
 - pydot>=1.2.3
 - cython
-- nipype>=1.1.6
+- nipype>=1.1.7
 
 - pip:
     - sphinx-argparse

fmriprep/__about__.py

@@ -88,7 +88,7 @@
     'indexed_gzip>=0.8.8',
     'nibabel>=2.2.1',
     'nilearn',
-    'nipype>=1.1.6',
+    'nipype>=1.1.7',
     'nitime',
     'niworkflows>=0.5.2.post5,<0.5.3',
     'numpy',

fmriprep/interfaces/patches.py

@@ -12,49 +12,104 @@
 
 from numpy.linalg.linalg import LinAlgError
 from nipype.algorithms import confounds as nac
+from nipype.interfaces.base import File
+from nipype.interfaces.mixins import reporting
 
 
-class RobustACompCor(nac.ACompCor):
-    """
-    Runs aCompCor several times if it suddenly fails with
-    https://github.com/poldracklab/fmriprep/issues/776
+class RetryCompCorInputSpecMixin(reporting.ReportCapableInputSpec):
+    out_report = File('report.html', usedefault=True, hash_files=False,
+                      desc='filename for warning HTML snippet')
 
-    """
 
+class RetryCompCorMixin(reporting.ReportCapableInterface):
     def _run_interface(self, runtime):
+        warn = self.inputs.failure_mode == 'NaN'
+
         failures = 0
+        save_exc = None
         while True:
+            success = True
+            # Identifiy success/failure in both error and NaN mode
             try:
-                runtime = super(RobustACompCor, self)._run_interface(runtime)
+                runtime = super()._run_interface(runtime)
+                if warn and self._is_allnans():
+                    success = False
+            except LinAlgError as exc:
+                success = False
+                save_exc = exc
+
+            if success:
                 break
-            except LinAlgError:
-                failures += 1
-                if failures > 10:
-                    raise
-                start = (failures - 1) * 10
-                sleep(randint(start + 4, start + 10))
+
+            failures += 1
+            if failures > 10:
+                if warn:
+                    break
+                raise save_exc
+            start = (failures - 1) * 10
+            sleep(randint(start + 4, start + 10))
 
         return runtime
 
+    def _is_allnans(self):
+        import numpy as np
+        outputs = self._list_outputs()
+        components = np.loadtxt(outputs['components_file'], skiprows=1)
+        return np.isnan(components).all()
+
+    def _generate_report(self):
+        snippet = '<!-- {} completed without error -->'.format(self._header)
+        if self._is_allnans():
+            snippet = '''\
+<p class="elem-desc">
+    Warning: {} components could not be estimated, due to a linear algebra error.
+    While not definitive, this may be an indication of a poor mask.
+    Please inspect the {} contours above to ensure that they are located
+    in the white matter/CSF.
+</p>
+'''.format(self._header, 'magenta' if self._header[0] == 'a' else 'blue')
+
+        with open(self._out_report, 'w') as fobj:
+            fobj.write(snippet)
 
-class RobustTCompCor(nac.TCompCor):
+
+class RobustACompCorInputSpec(RetryCompCorInputSpecMixin, nac.CompCorInputSpec):
+    pass
+
+
+class RobustACompCorOutputSpec(reporting.ReportCapableOutputSpec, nac.CompCorOutputSpec):
+    pass
+
+
+class RobustACompCor(RetryCompCorMixin, nac.ACompCor):
     """
-    Runs tCompCor several times if it suddenly fails with
-    https://github.com/poldracklab/fmriprep/issues/940
+    Runs aCompCor several times if it suddenly fails with
+    https://github.com/poldracklab/fmriprep/issues/776
+
+    Warns by default, rather than failing, on linear algebra errors.
+    https://github.com/poldracklab/fmriprep/issues/1433
 
     """
+    input_spec = RobustACompCorInputSpec
+    output_spec = RobustACompCorOutputSpec
 
-    def _run_interface(self, runtime):
-        failures = 0
-        while True:
-            try:
-                runtime = super(RobustTCompCor, self)._run_interface(runtime)
-                break
-            except LinAlgError:
-                failures += 1
-                if failures > 10:
-                    raise
-                start = (failures - 1) * 10
-                sleep(randint(start + 4, start + 10))
 
-        return runtime
+class RobustTCompCorInputSpec(RetryCompCorInputSpecMixin, nac.TCompCorInputSpec):
+    pass
+
+
+class RobustTCompCorOutputSpec(reporting.ReportCapableOutputSpec, nac.TCompCorOutputSpec):
+    pass
+
+
+class RobustTCompCor(RetryCompCorMixin, nac.TCompCor):
+    """
+    Runs tCompCor several times if it suddenly fails with
+    https://github.com/poldracklab/fmriprep/issues/776
+
+    Warns by default, rather than failing, on linear algebra errors.
+    https://github.com/poldracklab/fmriprep/issues/1433
+
+    """
+    input_spec = RobustTCompCorInputSpec
+    output_spec = RobustTCompCorOutputSpec

fmriprep/viz/config.json

@@ -98,6 +98,16 @@
                 "title": "ROIs in BOLD space",
                 "description": "Brain mask calculated on the BOLD signal (red contour), along with the masks used for a/tCompCor.<br />The aCompCor mask (magenta contour) is a conservative CSF and white-matter mask for extracting physiological and movement confounds. <br />The fCompCor mask (blue contour) contains the top 5% most variable voxels within a heavily-eroded brain-mask."
             },
+            {
+                "name": "epi/warn_tcompcor",
+                "file_pattern": "func/.*_acompcor\\.",
+                "raw": true
+            },
+            {
+                "name": "epi/warn_tcompcor",
+                "file_pattern": "func/.*_tcompcor\\.",
+                "raw": true
+            },
             {
                 "name": "epi_mean_t1_registration/flirt",
                 "file_pattern": "func/.*_flirtnobbr",

fmriprep/workflows/bold/confounds.py

@@ -183,12 +183,13 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"):
     # a/t-CompCor
     tcompcor = pe.Node(
         TCompCor(components_file='tcompcor.tsv', header_prefix='t_comp_cor_', pre_filter='cosine',
-                 save_pre_filter=True, percentile_threshold=.05),
+                 save_pre_filter=True, percentile_threshold=.05, failure_mode='NaN',
+                 generate_report=True),
         name="tcompcor", mem_gb=mem_gb)
 
     acompcor = pe.Node(
         ACompCor(components_file='acompcor.tsv', header_prefix='a_comp_cor_', pre_filter='cosine',
-                 save_pre_filter=True),
+                 save_pre_filter=True, failure_mode='NaN', generate_report=True),
         name="acompcor", mem_gb=mem_gb)
 
     # Set TR if present
@@ -223,6 +224,16 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"):
         name='ds_report_bold_rois', run_without_submitting=True,
         mem_gb=DEFAULT_MEMORY_MIN_GB)
 
+    ds_report_warn_acompcor = pe.Node(
+        DerivativesDataSink(suffix='acompcor'),
+        name='ds_report_warn_acompcor', run_without_submitting=True,
+        mem_gb=DEFAULT_MEMORY_MIN_GB)
+
+    ds_report_warn_tcompcor = pe.Node(
+        DerivativesDataSink(suffix='tcompcor'),
+        name='ds_report_warn_tcompcor', run_without_submitting=True,
+        mem_gb=DEFAULT_MEMORY_MIN_GB)
+
     def _pick_csf(files):
         return files[0]
 
@@ -303,6 +314,8 @@ def _pick_wm(files):
         (acc_msk, mrg_compcor, [('out', 'in2')]),
         (mrg_compcor, rois_plot, [('out', 'in_rois')]),
         (rois_plot, ds_report_bold_rois, [('out_report', 'in_file')]),
+        (acompcor, ds_report_warn_acompcor, [('out_report', 'in_file')]),
+        (tcompcor, ds_report_warn_tcompcor, [('out_report', 'in_file')]),
     ])
 
     return workflow