From d45dc87ac09891e80422ad44bbc98d609b3b2abd Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Thu, 15 Sep 2016 02:35:42 +0000
Subject: [PATCH 01/29] extractsignal inputs/outputs

---
 nipype/algorithms/stats.py | 79 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 79 insertions(+)
 create mode 100644 nipype/algorithms/stats.py

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
new file mode 100644
index 0000000000..1da1fb2fab
--- /dev/null
+++ b/nipype/algorithms/stats.py
@@ -0,0 +1,79 @@
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+'''
+Algorithms to compute statistics on :abbr:`fMRI (functional MRI)`
+
+    Change directory to provide relative paths for doctests
+    >>> import os
+    >>> filepath = os.path.dirname(os.path.realpath(__file__))
+    >>> datadir = os.path.realpath(os.path.join(filepath, '../testing/data'))
+    >>> os.chdir(datadir)
+
+'''
+from __future__ import print_function, division, unicode_literals, absolute_import
+from builtins import str, zip, range, open
+
+import os
+import os.path as op
+
+import nibabel as nb
+import numpy as np
+from scipy import linalg
+from scipy.special import legendre
+
+from .. import logging
+from ..external.due import due, Doi, BibTeX
+from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
+                               BaseInterfaceInputSpec, File, isdefined,
+                               InputMultiPath, ListStr)
+IFLOG = logging.getLogger('interface')
+
+class StatExtractionInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='4-D fMRI nii file')
+    label_file = File(exists=True, mandatory=False,
+                      desc='a 3-D label image, with 0 denoting background, or '
+                      'a 4-D file of probability maps. If this is not '
+                      'provided, this interface outputs one stat.')
+    out_file = File('stats.tsv', usedefault=True, exists=False, mandatory=False,
+                    desc='The name of the file to output the stats to. '
+                    'stats.tsv by default')
+    class_labels = ListStr(mandatory=False,
+                           desc='Human-readable labels for each segment in the '
+                           'label file, in order. The length of class_labels '
+                           'must be equal to or less than the number of '
+                           'segments.')
+    stat = Enum(('mean',), mandatory=False, default='mean', usedefault=True,
+                desc='The stat you wish to calculate on each segment. '
+                'The default is findig the mean')
+
+class StatExtractionOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='tsv file containing the computed stats, '
+                    'with as many stats as there are labels and as many rows '
+                    'as there are timepoints in in_file')
+
+class StatExtraction(BaseInterface):
+    '''
+    Extracts time series stats over tissue classes or brain regions
+
+    >>> seinterface = StatExtraction()
+    >>> seinterface.inputs.in_file = 'functional.nii'
+    >>> seinterface.inputs.in_file = 'segmentation0.nii'
+    >>> seinterface.inputs.out_file = 'means.tsv'
+    >>> segments =['background', 'CSF', 'gray', 'white']
+    >>> seinterface.inputs.class_labels = segments
+    >>> seinterface.inputs.stat = 'mean'
+    '''
+    input_spec = StatExtractionInputSpec
+    output_spec = StatExtractionOutputSpec
+
+    def _run_interface(self, runtime):
+        # assert/check inputs make sense
+
+        
+        pass
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs['out_file'] = self.inputs.out_file
+        return outputs

From b736086e519f53a0944a2942ada592cb42f2511e Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 01:29:39 +0000
Subject: [PATCH 02/29] complete signal extraction interface

---
 nipype/algorithms/stats.py              | 46 ++++++++++-----
 nipype/algorithms/tests/test_compcor.py | 13 ++---
 nipype/algorithms/tests/test_stats.py   | 74 +++++++++++++++++++++++++
 3 files changed, 111 insertions(+), 22 deletions(-)
 create mode 100644 nipype/algorithms/tests/test_stats.py

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 1da1fb2fab..6b9d5c7c3b 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -21,31 +21,33 @@
 import numpy as np
 from scipy import linalg
 from scipy.special import legendre
+import nilearn.input_data as nl
 
 from .. import logging
 from ..external.due import due, Doi, BibTeX
 from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
                                BaseInterfaceInputSpec, File, isdefined,
-                               InputMultiPath, ListStr)
+                               InputMultiPath)
 IFLOG = logging.getLogger('interface')
 
 class StatExtractionInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='4-D fMRI nii file')
-    label_file = File(exists=True, mandatory=False,
+    label_file = File(exists=True, mandatory=True,
                       desc='a 3-D label image, with 0 denoting background, or '
                       'a 4-D file of probability maps. If this is not '
                       'provided, this interface outputs one stat.')
+    class_labels = traits.List(mandatory=True,
+                               desc='Human-readable labels for each segment '
+                               'in the label file, in order. The length of '
+                               'class_labels must be equal to the number of '
+                               'segments (background excluded)')
     out_file = File('stats.tsv', usedefault=True, exists=False, mandatory=False,
                     desc='The name of the file to output the stats to. '
                     'stats.tsv by default')
-    class_labels = ListStr(mandatory=False,
-                           desc='Human-readable labels for each segment in the '
-                           'label file, in order. The length of class_labels '
-                           'must be equal to or less than the number of '
-                           'segments.')
-    stat = Enum(('mean',), mandatory=False, default='mean', usedefault=True,
-                desc='The stat you wish to calculate on each segment. '
-                'The default is findig the mean')
+    stat = traits.Enum(('mean',), mandatory=False, default='mean',
+                       usedefault=True,
+                       desc='The stat you wish to calculate on each segment. '
+                       'The default is finding the mean')
 
 class StatExtractionOutputSpec(TraitedSpec):
     out_file = File(exists=True, desc='tsv file containing the computed stats, '
@@ -58,9 +60,9 @@ class StatExtraction(BaseInterface):
 
     >>> seinterface = StatExtraction()
     >>> seinterface.inputs.in_file = 'functional.nii'
-    >>> seinterface.inputs.in_file = 'segmentation0.nii'
+    >>> seinterface.inputs.in_file = 'segmentation0.nii.gz'
     >>> seinterface.inputs.out_file = 'means.tsv'
-    >>> segments =['background', 'CSF', 'gray', 'white']
+    >>> segments = ['CSF', 'gray', 'white']
     >>> seinterface.inputs.class_labels = segments
     >>> seinterface.inputs.stat = 'mean'
     '''
@@ -68,10 +70,24 @@ class StatExtraction(BaseInterface):
     output_spec = StatExtractionOutputSpec
 
     def _run_interface(self, runtime):
-        # assert/check inputs make sense
+        ins = self.inputs
+        if ins.stat == 'mean': # always true for now
+            nlmasker = nl.NiftiLabelsMasker(ins.label_file).fit()
+            region_signals = nlmasker.transform_single_imgs(ins.in_file)
 
-        
-        pass
+            num_labels_found = region_signals.shape[1]
+            if len(ins.class_labels) != num_labels_found:
+                raise ValueError('The length of class_labels {} does not '
+                                 'match the number of regions {} found in '
+                                 'label_file {}'.format(ins.class_labels,
+                                                        num_labels_found,
+                                                        ins.label_file))
+
+            output = np.vstack((ins.class_labels, region_signals.astype(str)))
+
+            # save output
+            np.savetxt(ins.out_file, output, fmt=b'%s', delimiter='\t')
+        return runtime
 
     def _list_outputs(self):
         outputs = self._outputs().get()
diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py
index d7b81a8cd2..24e714b92e 100644
--- a/nipype/algorithms/tests/test_compcor.py
+++ b/nipype/algorithms/tests/test_compcor.py
@@ -14,19 +14,18 @@
 class TestCompCor(unittest.TestCase):
     ''' Note: Tests currently do a poor job of testing functionality '''
 
-    filenames = {'functionalnii': 'compcorfunc.nii',
-                 'masknii': 'compcormask.nii',
-                 'components_file': None}
+    filenames = {
+        'in_file': None,
+        'label_file': None,
+        'out_file': 'stats.tsv'
+    }
 
     def setUp(self):
         # setup
         self.orig_dir = os.getcwd()
         self.temp_dir = tempfile.mkdtemp()
         os.chdir(self.temp_dir)
-        noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape)
-        self.realigned_file = utils.save_toy_nii(self.fake_data + noise,
-                                                 self.filenames['functionalnii'])
-        mask = np.ones(self.fake_data.shape[:3])
+
         mask[0,0,0] = 0
         mask[0,0,1] = 0
         self.mask_file = utils.save_toy_nii(mask, self.filenames['masknii'])
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
new file mode 100644
index 0000000000..110cbf1250
--- /dev/null
+++ b/nipype/algorithms/tests/test_stats.py
@@ -0,0 +1,74 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+import nipype
+from ...testing import (assert_equal, assert_true, assert_false, skipif, utils,
+                        assert_almost_equal)
+from .. import stats
+
+import unittest
+import nibabel as nb
+import numpy as np
+import os
+import tempfile
+import shutil
+
+class TestStatExtraction(unittest.TestCase):
+
+    filenames = {
+        'in_file': 'fmri.nii',
+        'label_file': 'labels.nii',
+        'out_file': 'stats.tsv'
+    }
+
+    def setUp(self):
+        self.orig_dir = os.getcwd()
+        self.temp_dir = tempfile.mkdtemp()
+        os.chdir(self.temp_dir)
+
+        utils.save_toy_nii(self.fake_fmri_data, self.filenames['in_file'])
+        utils.save_toy_nii(self.fake_label_data, self.filenames['label_file'])
+
+    @skipif(True)
+    def test_stat_extraction(self):
+        # setup
+        wanted = [[]]
+        num_timepoints_wanted = self.fake_fmri_data.shape[3]
+        # run
+
+        labels_wanted = ['csf', 'gray', 'white']
+        stats.StatExtraction(in_file=self.filenames['in_file'],
+                             label_file=self.filenames['label_file'],
+                             class_labels=labels_wanted).run()
+        # assert
+        with open(self.filenames['out_file'], 'r') as output:
+            got = [line.split() for line in output]
+            labels_got = got.pop(0)
+            assert_equal(labels_got, labels_wanted)
+            assert_equal(len(got), num_timepoints_wanted)
+            for time in range(len(got)):
+                for segment in range(len(got[time])):
+                    assert_almost_equal(got[time][segment],
+                                        wanted[time][segment])
+
+    def tearDown(self):
+        os.chdir(self.orig_dir)
+        shutil.rmtree(self.temp_dir)
+
+    fake_fmri_data = np.array([[[[ 2, -1,  4, -2,  3],
+                                 [ 4, -2, -5, -1,  0]],
+                                
+                                [[-2,  0,  1,  4,  4],
+                                 [-5,  3, -3,  1, -5]]],
+                               
+                               
+                               [[[ 2, -2, -1, -2, -5],
+                                 [ 3,  0,  3, -5, -2]],
+                                
+                                [[-4, -2, -2,  1, -2],
+                                 [ 3,  1,  4, -3, -2]]]])
+    
+    fake_label_data = np.array([[[1, 0],
+                                 [3, 1]],
+
+                                [[2, 0],
+                                 [1, 3]]])

From e46bd0051b6ec2e0fb5ab510d02e363094690e3c Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 01:38:17 +0000
Subject: [PATCH 03/29] change name of interface

---
 nipype/algorithms/stats.py            | 30 +++++++++++++--------------
 nipype/algorithms/tests/test_stats.py | 12 +++++------
 2 files changed, 21 insertions(+), 21 deletions(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 6b9d5c7c3b..eed6071c2e 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -30,35 +30,35 @@
                                InputMultiPath)
 IFLOG = logging.getLogger('interface')
 
-class StatExtractionInputSpec(BaseInterfaceInputSpec):
+class SignalExtractionInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='4-D fMRI nii file')
     label_file = File(exists=True, mandatory=True,
                       desc='a 3-D label image, with 0 denoting background, or '
-                      'a 4-D file of probability maps. If this is not '
-                      'provided, this interface outputs one stat.')
+                      'a 4-D file of probability maps.')
     class_labels = traits.List(mandatory=True,
                                desc='Human-readable labels for each segment '
                                'in the label file, in order. The length of '
                                'class_labels must be equal to the number of '
                                'segments (background excluded)')
-    out_file = File('stats.tsv', usedefault=True, exists=False, mandatory=False,
-                    desc='The name of the file to output the stats to. '
-                    'stats.tsv by default')
+    out_file = File('signals.tsv', usedefault=True, exists=False,
+                    mandatory=False, desc='The name of the file to output to. '
+                    'signals.tsv by default')
     stat = traits.Enum(('mean',), mandatory=False, default='mean',
                        usedefault=True,
                        desc='The stat you wish to calculate on each segment. '
                        'The default is finding the mean')
 
-class StatExtractionOutputSpec(TraitedSpec):
-    out_file = File(exists=True, desc='tsv file containing the computed stats, '
-                    'with as many stats as there are labels and as many rows '
-                    'as there are timepoints in in_file')
+class SignalExtractionOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='tsv file containing the computed '
+                    'signals, with as many columns as there are labels and as '
+                    'many rows as there are timepoints in in_file, plus a '
+                    'header row with values from class_labels')
 
-class StatExtraction(BaseInterface):
+class SignalExtraction(BaseInterface):
     '''
-    Extracts time series stats over tissue classes or brain regions
+    Extracts signals over tissue classes or brain regions
 
-    >>> seinterface = StatExtraction()
+    >>> seinterface = SignalExtraction()
     >>> seinterface.inputs.in_file = 'functional.nii'
     >>> seinterface.inputs.in_file = 'segmentation0.nii.gz'
     >>> seinterface.inputs.out_file = 'means.tsv'
@@ -66,8 +66,8 @@ class StatExtraction(BaseInterface):
     >>> seinterface.inputs.class_labels = segments
     >>> seinterface.inputs.stat = 'mean'
     '''
-    input_spec = StatExtractionInputSpec
-    output_spec = StatExtractionOutputSpec
+    input_spec = SignalExtractionInputSpec
+    output_spec = SignalExtractionOutputSpec
 
     def _run_interface(self, runtime):
         ins = self.inputs
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index 110cbf1250..6cd1213369 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -12,12 +12,12 @@
 import tempfile
 import shutil
 
-class TestStatExtraction(unittest.TestCase):
+class TestSignalExtraction(unittest.TestCase):
 
     filenames = {
         'in_file': 'fmri.nii',
         'label_file': 'labels.nii',
-        'out_file': 'stats.tsv'
+        'out_file': 'signals.tsv'
     }
 
     def setUp(self):
@@ -29,16 +29,16 @@ def setUp(self):
         utils.save_toy_nii(self.fake_label_data, self.filenames['label_file'])
 
     @skipif(True)
-    def test_stat_extraction(self):
+    def test_signal_extraction(self):
         # setup
         wanted = [[]]
         num_timepoints_wanted = self.fake_fmri_data.shape[3]
         # run
 
         labels_wanted = ['csf', 'gray', 'white']
-        stats.StatExtraction(in_file=self.filenames['in_file'],
-                             label_file=self.filenames['label_file'],
-                             class_labels=labels_wanted).run()
+        stats.SignalExtraction(in_file=self.filenames['in_file'],
+                               label_file=self.filenames['label_file'],
+                               class_labels=labels_wanted).run()
         # assert
         with open(self.filenames['out_file'], 'r') as output:
             got = [line.split() for line in output]

From d51ab4326cf0e22cba32c28928fab5d7df86770f Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 01:50:31 +0000
Subject: [PATCH 04/29] test failure case

---
 nipype/algorithms/tests/test_stats.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index 6cd1213369..d295f6a2a9 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -2,7 +2,7 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 import nipype
 from ...testing import (assert_equal, assert_true, assert_false, skipif, utils,
-                        assert_almost_equal)
+                        assert_almost_equal, raises)
 from .. import stats
 
 import unittest
@@ -50,6 +50,13 @@ def test_signal_extraction(self):
                     assert_almost_equal(got[time][segment],
                                         wanted[time][segment])
 
+    @raises(ValueError)
+    def test_signal_extraction_bad_class_labels(self):
+        # run
+        stats.SignalExtraction(in_file=self.filenames['in_file'],
+                               label_file=self.filenames['label_file'],
+                               class_labels=['bad']).run()
+
     def tearDown(self):
         os.chdir(self.orig_dir)
         shutil.rmtree(self.temp_dir)

From 02faccc1898cfe63d8562508bbfa392dccec3540 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 02:13:16 +0000
Subject: [PATCH 05/29] finalize test

---
 nipype/algorithms/tests/test_stats.py | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index d295f6a2a9..ff6e38b6e6 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -28,10 +28,13 @@ def setUp(self):
         utils.save_toy_nii(self.fake_fmri_data, self.filenames['in_file'])
         utils.save_toy_nii(self.fake_label_data, self.filenames['label_file'])
 
-    @skipif(True)
     def test_signal_extraction(self):
         # setup
-        wanted = [[]]
+        wanted = [[-2.33333, 2, .5],
+                  [0, -2, .5],
+                  [-.3333333, -1, 2.5],
+                  [0, -2, .5],
+                  [-1.3333333, -5, 1]]
         num_timepoints_wanted = self.fake_fmri_data.shape[3]
         # run
 
@@ -42,13 +45,16 @@ def test_signal_extraction(self):
         # assert
         with open(self.filenames['out_file'], 'r') as output:
             got = [line.split() for line in output]
-            labels_got = got.pop(0)
+            labels_got = got.pop(0) # remove header
             assert_equal(labels_got, labels_wanted)
             assert_equal(len(got), num_timepoints_wanted)
+            # convert from string to float
+            got = [[float(num) for num in row] for row in got]
             for time in range(len(got)):
+                assert_equal(len(labels_wanted), len(got[time]))
                 for segment in range(len(got[time])):
                     assert_almost_equal(got[time][segment],
-                                        wanted[time][segment])
+                                        wanted[time][segment], decimal=1)
 
     @raises(ValueError)
     def test_signal_extraction_bad_class_labels(self):

From 5709f055a3c53a56bcc54f144cd5a706b9a95fa6 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 02:17:05 +0000
Subject: [PATCH 06/29] pep8

---
 nipype/algorithms/tests/test_stats.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index ff6e38b6e6..e92603120a 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -69,17 +69,17 @@ def tearDown(self):
 
     fake_fmri_data = np.array([[[[ 2, -1,  4, -2,  3],
                                  [ 4, -2, -5, -1,  0]],
-                                
+
                                 [[-2,  0,  1,  4,  4],
                                  [-5,  3, -3,  1, -5]]],
-                               
-                               
+
+
                                [[[ 2, -2, -1, -2, -5],
                                  [ 3,  0,  3, -5, -2]],
-                                
+
                                 [[-4, -2, -2,  1, -2],
                                  [ 3,  1,  4, -3, -2]]]])
-    
+
     fake_label_data = np.array([[[1, 0],
                                  [3, 1]],
 

From 0a5933308993c921f8251132cddc0b4381e19785 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 02:37:51 +0000
Subject: [PATCH 07/29] unused imports, pep8

---
 nipype/algorithms/stats.py            | 15 +++--------
 nipype/algorithms/tests/test_stats.py | 37 ++++++++++++---------------
 2 files changed, 21 insertions(+), 31 deletions(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index eed6071c2e..9d5c09b6e7 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -11,23 +11,16 @@
     >>> os.chdir(datadir)
 
 '''
-from __future__ import print_function, division, unicode_literals, absolute_import
-from builtins import str, zip, range, open
+from __future__ import (print_function, division, unicode_literals,
+                        absolute_import)
+from builtins import str
 
-import os
-import os.path as op
-
-import nibabel as nb
 import numpy as np
-from scipy import linalg
-from scipy.special import legendre
 import nilearn.input_data as nl
 
 from .. import logging
-from ..external.due import due, Doi, BibTeX
 from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
-                               BaseInterfaceInputSpec, File, isdefined,
-                               InputMultiPath)
+                               BaseInterfaceInputSpec, File)
 IFLOG = logging.getLogger('interface')
 
 class SignalExtractionInputSpec(BaseInterfaceInputSpec):
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index e92603120a..3e91c2d48f 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -1,17 +1,15 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
-import nipype
-from ...testing import (assert_equal, assert_true, assert_false, skipif, utils,
-                        assert_almost_equal, raises)
-from .. import stats
-
 import unittest
-import nibabel as nb
-import numpy as np
 import os
 import tempfile
 import shutil
 
+import numpy as np
+
+from ...testing import (assert_equal, utils, assert_almost_equal, raises)
+from .. import stats
+
 class TestSignalExtraction(unittest.TestCase):
 
     filenames = {
@@ -50,11 +48,10 @@ def test_signal_extraction(self):
             assert_equal(len(got), num_timepoints_wanted)
             # convert from string to float
             got = [[float(num) for num in row] for row in got]
-            for time in range(len(got)):
-                assert_equal(len(labels_wanted), len(got[time]))
-                for segment in range(len(got[time])):
-                    assert_almost_equal(got[time][segment],
-                                        wanted[time][segment], decimal=1)
+            for i, time in enumerate(got):
+                assert_equal(len(labels_wanted), len(time))
+                for j, segment in enumerate(time):
+                    assert_almost_equal(segment, wanted[i][j], decimal=1)
 
     @raises(ValueError)
     def test_signal_extraction_bad_class_labels(self):
@@ -67,18 +64,18 @@ def tearDown(self):
         os.chdir(self.orig_dir)
         shutil.rmtree(self.temp_dir)
 
-    fake_fmri_data = np.array([[[[ 2, -1,  4, -2,  3],
-                                 [ 4, -2, -5, -1,  0]],
+    fake_fmri_data = np.array([[[[2, -1, 4, -2, 3],
+                                 [4, -2, -5, -1, 0]],
 
-                                [[-2,  0,  1,  4,  4],
-                                 [-5,  3, -3,  1, -5]]],
+                                [[-2, 0, 1, 4, 4],
+                                 [-5, 3, -3, 1, -5]]],
 
 
-                               [[[ 2, -2, -1, -2, -5],
-                                 [ 3,  0,  3, -5, -2]],
+                               [[[2, -2, -1, -2, -5],
+                                 [3, 0, 3, -5, -2]],
 
-                                [[-4, -2, -2,  1, -2],
-                                 [ 3,  1,  4, -3, -2]]]])
+                                [[-4, -2, -2, 1, -2],
+                                 [3, 1, 4, -3, -2]]]])
 
     fake_label_data = np.array([[[1, 0],
                                  [3, 1]],

From e0faaf5c7587a2bb25372b696ab9bb74ce7bf726 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 02:42:37 +0000
Subject: [PATCH 08/29] undo accidental changes to test_compcor

---
 nipype/algorithms/tests/test_compcor.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py
index 24e714b92e..d7b81a8cd2 100644
--- a/nipype/algorithms/tests/test_compcor.py
+++ b/nipype/algorithms/tests/test_compcor.py
@@ -14,18 +14,19 @@
 class TestCompCor(unittest.TestCase):
     ''' Note: Tests currently do a poor job of testing functionality '''
 
-    filenames = {
-        'in_file': None,
-        'label_file': None,
-        'out_file': 'stats.tsv'
-    }
+    filenames = {'functionalnii': 'compcorfunc.nii',
+                 'masknii': 'compcormask.nii',
+                 'components_file': None}
 
     def setUp(self):
         # setup
         self.orig_dir = os.getcwd()
         self.temp_dir = tempfile.mkdtemp()
         os.chdir(self.temp_dir)
-
+        noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape)
+        self.realigned_file = utils.save_toy_nii(self.fake_data + noise,
+                                                 self.filenames['functionalnii'])
+        mask = np.ones(self.fake_data.shape[:3])
         mask[0,0,0] = 0
         mask[0,0,1] = 0
         self.mask_file = utils.save_toy_nii(mask, self.filenames['masknii'])

From 160994694deac884e25c41fe346e9e4881dc706a Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 02:56:43 +0000
Subject: [PATCH 09/29] add nilearn to dependencies

---
 requirements.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/requirements.txt b/requirements.txt
index ef66036744..69bf80f8a3 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -12,3 +12,4 @@ xvfbwrapper
 psutil
 funcsigs
 configparser
+nilearn

From fffd9b56737dc4653860344ac34f432edbab0db3 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 19:58:58 +0000
Subject: [PATCH 10/29] add sklearn explicitly to dependencies like
 https://github.com/poldracklab/fmriprep/commit/92830add5ad54f34325a785d9479abf0c9030d96

---
 requirements.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/requirements.txt b/requirements.txt
index 69bf80f8a3..b8aa8aed2f 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -13,3 +13,4 @@ psutil
 funcsigs
 configparser
 nilearn
+sklearn

From 700ecce40d4694c9557d68d805a573e3b884b77e Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 21:08:06 +0000
Subject: [PATCH 11/29] add detrend param, doc

---
 nipype/algorithms/stats.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 9d5c09b6e7..7e2fcce6d2 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -32,7 +32,9 @@ class SignalExtractionInputSpec(BaseInterfaceInputSpec):
                                desc='Human-readable labels for each segment '
                                'in the label file, in order. The length of '
                                'class_labels must be equal to the number of '
-                               'segments (background excluded)')
+                               'segments (background excluded). This list '
+                               'corresponds to the class labels in label_file '
+                               'in ascending order')
     out_file = File('signals.tsv', usedefault=True, exists=False,
                     mandatory=False, desc='The name of the file to output to. '
                     'signals.tsv by default')
@@ -40,6 +42,8 @@ class SignalExtractionInputSpec(BaseInterfaceInputSpec):
                        usedefault=True,
                        desc='The stat you wish to calculate on each segment. '
                        'The default is finding the mean')
+    detrend = traits.Bool(False, usedefault=True, mandatory=False,
+                          desc='If True, perform detrending using nilearn.')
 
 class SignalExtractionOutputSpec(TraitedSpec):
     out_file = File(exists=True, desc='tsv file containing the computed '
@@ -64,8 +68,11 @@ class SignalExtraction(BaseInterface):
 
     def _run_interface(self, runtime):
         ins = self.inputs
+
         if ins.stat == 'mean': # always true for now
-            nlmasker = nl.NiftiLabelsMasker(ins.label_file).fit()
+            nlmasker = nl.NiftiLabelsMasker(ins.label_file,
+                                            detrend=ins.detrend)
+            nlmasker.fit()
             region_signals = nlmasker.transform_single_imgs(ins.in_file)
 
             num_labels_found = region_signals.shape[1]

From 4f79c18911ec26acc1b914e18be0d377562cef7a Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Fri, 16 Sep 2016 21:21:18 +0000
Subject: [PATCH 12/29] make nilearn optional

---
 nipype/algorithms/stats.py            |  3 ++-
 nipype/algorithms/tests/test_stats.py | 12 +++++++++++-
 requirements.txt                      |  2 --
 3 files changed, 13 insertions(+), 4 deletions(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 7e2fcce6d2..9c8187276d 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -16,7 +16,6 @@
 from builtins import str
 
 import numpy as np
-import nilearn.input_data as nl
 
 from .. import logging
 from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
@@ -67,6 +66,8 @@ class SignalExtraction(BaseInterface):
     output_spec = SignalExtractionOutputSpec
 
     def _run_interface(self, runtime):
+        import nilearn.input_data as nl
+
         ins = self.inputs
 
         if ins.stat == 'mean': # always true for now
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index 3e91c2d48f..c1479d7501 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -7,9 +7,17 @@
 
 import numpy as np
 
-from ...testing import (assert_equal, utils, assert_almost_equal, raises)
+from ...testing import (assert_equal, utils, assert_almost_equal, raises,
+                        skipif)
 from .. import stats
 
+no_nilearn = True
+try:
+    import nilearn
+    no_nilearn = False
+except ImportError:
+    pass
+
 class TestSignalExtraction(unittest.TestCase):
 
     filenames = {
@@ -26,6 +34,7 @@ def setUp(self):
         utils.save_toy_nii(self.fake_fmri_data, self.filenames['in_file'])
         utils.save_toy_nii(self.fake_label_data, self.filenames['label_file'])
 
+    @skipif(no_nilearn)
     def test_signal_extraction(self):
         # setup
         wanted = [[-2.33333, 2, .5],
@@ -53,6 +62,7 @@ def test_signal_extraction(self):
                 for j, segment in enumerate(time):
                     assert_almost_equal(segment, wanted[i][j], decimal=1)
 
+    @skipif(no_nilearn)
     @raises(ValueError)
     def test_signal_extraction_bad_class_labels(self):
         # run
diff --git a/requirements.txt b/requirements.txt
index b8aa8aed2f..ef66036744 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -12,5 +12,3 @@ xvfbwrapper
 psutil
 funcsigs
 configparser
-nilearn
-sklearn

From b11c5f5801d29a77285aa6f17a7ac8dbf4352cdc Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Sat, 17 Sep 2016 02:02:22 +0000
Subject: [PATCH 13/29] refactoring & testing in preparation for 4d label
 files; pep8

---
 nipype/algorithms/stats.py            | 14 ++++++++------
 nipype/algorithms/tests/test_stats.py | 21 ++++++++++++++++++++-
 2 files changed, 28 insertions(+), 7 deletions(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 9c8187276d..b6900cd031 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -13,7 +13,6 @@
 '''
 from __future__ import (print_function, division, unicode_literals,
                         absolute_import)
-from builtins import str
 
 import numpy as np
 
@@ -71,11 +70,7 @@ def _run_interface(self, runtime):
         ins = self.inputs
 
         if ins.stat == 'mean': # always true for now
-            nlmasker = nl.NiftiLabelsMasker(ins.label_file,
-                                            detrend=ins.detrend)
-            nlmasker.fit()
-            region_signals = nlmasker.transform_single_imgs(ins.in_file)
-
+            region_signals = self._3d_label_handler(nl, ins)
             num_labels_found = region_signals.shape[1]
             if len(ins.class_labels) != num_labels_found:
                 raise ValueError('The length of class_labels {} does not '
@@ -90,6 +85,13 @@ def _run_interface(self, runtime):
             np.savetxt(ins.out_file, output, fmt=b'%s', delimiter='\t')
         return runtime
 
+    def _3d_label_handler(self, nl, ins):
+        nlmasker = nl.NiftiLabelsMasker(ins.label_file,
+                                        detrend=ins.detrend)
+        nlmasker.fit()
+        region_signals = nlmasker.transform_single_imgs(ins.in_file)
+        return region_signals
+
     def _list_outputs(self):
         outputs = self._outputs().get()
         outputs['out_file'] = self.inputs.out_file
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index c1479d7501..081c97dff2 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -13,7 +13,7 @@
 
 no_nilearn = True
 try:
-    import nilearn
+    __import__('nilearn')
     no_nilearn = False
 except ImportError:
     pass
@@ -23,6 +23,7 @@ class TestSignalExtraction(unittest.TestCase):
     filenames = {
         'in_file': 'fmri.nii',
         'label_file': 'labels.nii',
+        '4d_label_file': '4dlabels.nii',
         'out_file': 'signals.tsv'
     }
 
@@ -70,6 +71,11 @@ def test_signal_extraction_bad_class_labels(self):
                                label_file=self.filenames['label_file'],
                                class_labels=['bad']).run()
 
+    @skipif(no_nilearn)
+    def test_signal_extraction_4d_label(self):
+        # setup
+        utils.save_toy_nii(self.fake_4d_label_data, self.filenames['4d_label_file'])
+
     def tearDown(self):
         os.chdir(self.orig_dir)
         shutil.rmtree(self.temp_dir)
@@ -92,3 +98,16 @@ def tearDown(self):
 
                                 [[2, 0],
                                  [1, 3]]])
+
+    fake_4d_label_data = np.array([[[[0.2, 0.3, 0.5],
+                                     [0.1, 0.1, 0.8]],
+
+                                    [[0.1, 0.3, 0.6],
+                                     [0.3, 0.4, 0.3]]],
+
+
+                                   [[[0.2, 0.2, 0.6],
+                                     [0., 0.3, 0.7]],
+
+                                    [[0.3, 0.3, 0.4],
+                                     [0.3, 0.4, 0.3]]]])

From fb0d6b9fa7f8306821ebe9528aec0c5d80aff352 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Sun, 18 Sep 2016 00:17:37 +0000
Subject: [PATCH 14/29] nilearn method

---
 nipype/algorithms/stats.py            | 20 +++++--
 nipype/algorithms/tests/test_stats.py | 78 ++++++++++++++++++---------
 2 files changed, 68 insertions(+), 30 deletions(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index b6900cd031..3826a0a109 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -15,6 +15,7 @@
                         absolute_import)
 
 import numpy as np
+import nibabel as nb
 
 from .. import logging
 from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
@@ -68,9 +69,13 @@ def _run_interface(self, runtime):
         import nilearn.input_data as nl
 
         ins = self.inputs
+        labels = nb.load(ins.label_file)
 
         if ins.stat == 'mean': # always true for now
-            region_signals = self._3d_label_handler(nl, ins)
+            if len(labels.get_data().shape) == 3:
+                region_signals = self._3d_label_handler(nl, labels)
+            else:
+                region_signals = self._4d_label_handler(nl, labels)
             num_labels_found = region_signals.shape[1]
             if len(ins.class_labels) != num_labels_found:
                 raise ValueError('The length of class_labels {} does not '
@@ -85,11 +90,16 @@ def _run_interface(self, runtime):
             np.savetxt(ins.out_file, output, fmt=b'%s', delimiter='\t')
         return runtime
 
-    def _3d_label_handler(self, nl, ins):
-        nlmasker = nl.NiftiLabelsMasker(ins.label_file,
-                                        detrend=ins.detrend)
+    def _3d_label_handler(self, nl, labels):
+        nlmasker = nl.NiftiLabelsMasker(labels, detrend=self.inputs.detrend)
         nlmasker.fit()
-        region_signals = nlmasker.transform_single_imgs(ins.in_file)
+        region_signals = nlmasker.transform_single_imgs(self.inputs.in_file)
+        return region_signals
+
+    def _4d_label_handler(self, nl, labels):
+        nmmasker = nl.NiftiMapsMasker(labels, detrend=self.inputs.detrend)
+        nmmasker.fit()
+        region_signals = nmmasker.transform_single_imgs(self.inputs.in_file)
         return region_signals
 
     def _list_outputs(self):
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index 081c97dff2..40537fa5c4 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -26,6 +26,7 @@ class TestSignalExtraction(unittest.TestCase):
         '4d_label_file': '4dlabels.nii',
         'out_file': 'signals.tsv'
     }
+    labels = ['csf', 'gray', 'white']
 
     def setUp(self):
         self.orig_dir = os.getcwd()
@@ -37,44 +38,57 @@ def setUp(self):
 
     @skipif(no_nilearn)
     def test_signal_extraction(self):
-        # setup
-        wanted = [[-2.33333, 2, .5],
-                  [0, -2, .5],
-                  [-.3333333, -1, 2.5],
-                  [0, -2, .5],
-                  [-1.3333333, -5, 1]]
-        num_timepoints_wanted = self.fake_fmri_data.shape[3]
         # run
-
-        labels_wanted = ['csf', 'gray', 'white']
         stats.SignalExtraction(in_file=self.filenames['in_file'],
                                label_file=self.filenames['label_file'],
-                               class_labels=labels_wanted).run()
+                               class_labels=self.labels).run()
         # assert
-        with open(self.filenames['out_file'], 'r') as output:
-            got = [line.split() for line in output]
-            labels_got = got.pop(0) # remove header
-            assert_equal(labels_got, labels_wanted)
-            assert_equal(len(got), num_timepoints_wanted)
-            # convert from string to float
-            got = [[float(num) for num in row] for row in got]
-            for i, time in enumerate(got):
-                assert_equal(len(labels_wanted), len(time))
-                for j, segment in enumerate(time):
-                    assert_almost_equal(segment, wanted[i][j], decimal=1)
+        self.assert_expected_output(self.base_wanted)
 
     @skipif(no_nilearn)
     @raises(ValueError)
-    def test_signal_extraction_bad_class_labels(self):
+    def test_signal_extraction_bad_label_list(self):
         # run
         stats.SignalExtraction(in_file=self.filenames['in_file'],
                                label_file=self.filenames['label_file'],
                                class_labels=['bad']).run()
 
     @skipif(no_nilearn)
-    def test_signal_extraction_4d_label(self):
+    def test_signal_extraction_equiv_4d(self):
+        self._test_4d_label(self.base_wanted, self.fake_equiv_4d_label_data)
+
+    def test_signal_extraction_4d_(self):
+        self._test_4d_label([[-5.0652173913, -5.44565217391, 5.50543478261],
+                             [0, -2, .5],
+                             [-.3333333, -1, 2.5],
+                             [0, -2, .5],
+                             [-1.3333333, -5, 1]], self.fake_4d_label_data)
+
+    def _test_4d_label(self, wanted, fake_labels):
         # setup
-        utils.save_toy_nii(self.fake_4d_label_data, self.filenames['4d_label_file'])
+        utils.save_toy_nii(fake_labels, self.filenames['4d_label_file'])
+
+        # run
+        stats.SignalExtraction(in_file=self.filenames['in_file'],
+                               label_file=self.filenames['4d_label_file'],
+                               class_labels=self.labels).run()
+
+        self.assert_expected_output(wanted)
+
+    def assert_expected_output(self, wanted):
+        with open(self.filenames['out_file'], 'r') as output:
+            got = [line.split() for line in output]
+            labels_got = got.pop(0) # remove header
+            assert_equal(labels_got, self.labels)
+            assert_equal(len(got), self.fake_fmri_data.shape[3],
+                         'num rows and num volumes')
+            # convert from string to float
+            got = [[float(num) for num in row] for row in got]
+            for i, time in enumerate(got):
+                assert_equal(len(self.labels), len(time))
+                for j, segment in enumerate(time):
+                    assert_almost_equal(segment, wanted[i][j], decimal=1)
+
 
     def tearDown(self):
         os.chdir(self.orig_dir)
@@ -99,13 +113,27 @@ def tearDown(self):
                                 [[2, 0],
                                  [1, 3]]])
 
+    fake_equiv_4d_label_data = np.array([[[[1., 0., 0.],
+                                           [0., 0., 0.]],
+                                          [[0., 0., 1.],
+                                           [1., 0., 0.]]],
+                                         [[[0., 1., 0.],
+                                           [0., 0., 0.]],
+                                          [[1., 0., 0.],
+                                           [0., 0., 1.]]]])
+
+    base_wanted = [[-2.33333, 2, .5],
+                   [0, -2, .5],
+                   [-.3333333, -1, 2.5],
+                   [0, -2, .5],
+                   [-1.3333333, -5, 1]]
+
     fake_4d_label_data = np.array([[[[0.2, 0.3, 0.5],
                                      [0.1, 0.1, 0.8]],
 
                                     [[0.1, 0.3, 0.6],
                                      [0.3, 0.4, 0.3]]],
 
-
                                    [[[0.2, 0.2, 0.6],
                                      [0., 0.3, 0.7]],
 

From c1b3a128b27ff5e8fd348e47ba86a350ebf6e66f Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Sun, 18 Sep 2016 05:55:46 +0000
Subject: [PATCH 15/29] m

---
 nipype/algorithms/tests/test_stats.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index 40537fa5c4..534743a892 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -57,6 +57,7 @@ def test_signal_extraction_bad_label_list(self):
     def test_signal_extraction_equiv_4d(self):
         self._test_4d_label(self.base_wanted, self.fake_equiv_4d_label_data)
 
+    @skipif(no_nilearn)
     def test_signal_extraction_4d_(self):
         self._test_4d_label([[-5.0652173913, -5.44565217391, 5.50543478261],
                              [0, -2, .5],

From bc48a7433d9174a2af06adac10949a1cf810484c Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Sun, 18 Sep 2016 22:20:37 +0000
Subject: [PATCH 16/29] refactor/allow lists of 3d label files

---
 nipype/algorithms/stats.py            | 90 ++++++++++++++-------------
 nipype/algorithms/tests/test_stats.py | 20 +++---
 2 files changed, 58 insertions(+), 52 deletions(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 3826a0a109..0c2ee78d05 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -19,14 +19,16 @@
 
 from .. import logging
 from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
-                               BaseInterfaceInputSpec, File)
+                               BaseInterfaceInputSpec, File, InputMultiPath)
 IFLOG = logging.getLogger('interface')
 
 class SignalExtractionInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='4-D fMRI nii file')
-    label_file = File(exists=True, mandatory=True,
-                      desc='a 3-D label image, with 0 denoting background, or '
-                      'a 4-D file of probability maps.')
+    label_files = InputMultiPath(File(exists=True), mandatory=True,
+                                desc='a 3-D label image, with 0 denoting '
+                                'background, or a list of 3-D probability '
+                                'maps (one per label) or the equivalent 4D '
+                                'file.')
     class_labels = traits.List(mandatory=True,
                                desc='Human-readable labels for each segment '
                                'in the label file, in order. The length of '
@@ -37,10 +39,6 @@ class SignalExtractionInputSpec(BaseInterfaceInputSpec):
     out_file = File('signals.tsv', usedefault=True, exists=False,
                     mandatory=False, desc='The name of the file to output to. '
                     'signals.tsv by default')
-    stat = traits.Enum(('mean',), mandatory=False, default='mean',
-                       usedefault=True,
-                       desc='The stat you wish to calculate on each segment. '
-                       'The default is finding the mean')
     detrend = traits.Bool(False, usedefault=True, mandatory=False,
                           desc='If True, perform detrending using nilearn.')
 
@@ -56,51 +54,59 @@ class SignalExtraction(BaseInterface):
 
     >>> seinterface = SignalExtraction()
     >>> seinterface.inputs.in_file = 'functional.nii'
-    >>> seinterface.inputs.in_file = 'segmentation0.nii.gz'
+    >>> seinterface.inputs.label_files = 'segmentation0.nii.gz'
     >>> seinterface.inputs.out_file = 'means.tsv'
     >>> segments = ['CSF', 'gray', 'white']
     >>> seinterface.inputs.class_labels = segments
-    >>> seinterface.inputs.stat = 'mean'
+    >>> seinterface.inputs.detrend = True
     '''
     input_spec = SignalExtractionInputSpec
     output_spec = SignalExtractionOutputSpec
 
     def _run_interface(self, runtime):
-        import nilearn.input_data as nl
+        masker = self._process_inputs()
+
+        region_signals = masker.fit_transform(self.inputs.in_file)
+
+        output = np.vstack((self.inputs.class_labels, region_signals.astype(str)))
 
-        ins = self.inputs
-        labels = nb.load(ins.label_file)
-
-        if ins.stat == 'mean': # always true for now
-            if len(labels.get_data().shape) == 3:
-                region_signals = self._3d_label_handler(nl, labels)
-            else:
-                region_signals = self._4d_label_handler(nl, labels)
-            num_labels_found = region_signals.shape[1]
-            if len(ins.class_labels) != num_labels_found:
-                raise ValueError('The length of class_labels {} does not '
-                                 'match the number of regions {} found in '
-                                 'label_file {}'.format(ins.class_labels,
-                                                        num_labels_found,
-                                                        ins.label_file))
-
-            output = np.vstack((ins.class_labels, region_signals.astype(str)))
-
-            # save output
-            np.savetxt(ins.out_file, output, fmt=b'%s', delimiter='\t')
+        # save output
+        np.savetxt(self.inputs.out_file, output, fmt=b'%s', delimiter='\t')
         return runtime
 
-    def _3d_label_handler(self, nl, labels):
-        nlmasker = nl.NiftiLabelsMasker(labels, detrend=self.inputs.detrend)
-        nlmasker.fit()
-        region_signals = nlmasker.transform_single_imgs(self.inputs.in_file)
-        return region_signals
-
-    def _4d_label_handler(self, nl, labels):
-        nmmasker = nl.NiftiMapsMasker(labels, detrend=self.inputs.detrend)
-        nmmasker.fit()
-        region_signals = nmmasker.transform_single_imgs(self.inputs.in_file)
-        return region_signals
+    def _process_inputs(self):
+        ''' validate and  process inputs into useful form '''
+
+        import nilearn.input_data as nl
+
+        # determine form of label files, choose appropriate nilearn masker
+        if len(self.inputs.label_files) > 1: # list of 3D nifti images
+            masker = nl.NiftiMapsMasker(self.inputs.label_files)
+            n_labels = len(self.inputs.label_files)
+        else: # list of size one, containing either a 3d or a 4d file
+            label_data = nb.load(self.inputs.label_files[0])
+            if len(label_data.shape) == 4: # 4d file
+                masker = nl.NiftiMapsMasker(label_data)
+                n_labels = label_data.shape[3]
+            else: # 3d file
+                if np.amax(label_data) > 1: # 3d label file
+                    masker = nl.NiftiLabelsMasker(label_data)
+                    # assuming consecutive positive integers for regions
+                    n_labels = np.amax(label_data.get_data())
+                else: # most probably a single probability map for one label
+                    masker = nl.NiftiMapsMasker(label_data)
+                    n_labels = 1
+
+        # check label list size
+        if len(self.inputs.class_labels) != n_labels:
+            raise ValueError('The length of class_labels {} does not '
+                             'match the number of regions {} found in '
+                             'label_files {}'.format(self.inputs.class_labels,
+                                                    n_labels,
+                                                    self.inputs.label_files))
+
+        masker.set_params(detrend=self.inputs.detrend)
+        return masker
 
     def _list_outputs(self):
         outputs = self._outputs().get()
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index 534743a892..e198d087d5 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -22,7 +22,7 @@ class TestSignalExtraction(unittest.TestCase):
 
     filenames = {
         'in_file': 'fmri.nii',
-        'label_file': 'labels.nii',
+        'label_files': 'labels.nii',
         '4d_label_file': '4dlabels.nii',
         'out_file': 'signals.tsv'
     }
@@ -34,13 +34,13 @@ def setUp(self):
         os.chdir(self.temp_dir)
 
         utils.save_toy_nii(self.fake_fmri_data, self.filenames['in_file'])
-        utils.save_toy_nii(self.fake_label_data, self.filenames['label_file'])
+        utils.save_toy_nii(self.fake_label_data, self.filenames['label_files'])
 
     @skipif(no_nilearn)
     def test_signal_extraction(self):
         # run
         stats.SignalExtraction(in_file=self.filenames['in_file'],
-                               label_file=self.filenames['label_file'],
+                               label_files=self.filenames['label_files'],
                                class_labels=self.labels).run()
         # assert
         self.assert_expected_output(self.base_wanted)
@@ -50,7 +50,7 @@ def test_signal_extraction(self):
     def test_signal_extraction_bad_label_list(self):
         # run
         stats.SignalExtraction(in_file=self.filenames['in_file'],
-                               label_file=self.filenames['label_file'],
+                               label_files=self.filenames['label_files'],
                                class_labels=['bad']).run()
 
     @skipif(no_nilearn)
@@ -58,12 +58,12 @@ def test_signal_extraction_equiv_4d(self):
         self._test_4d_label(self.base_wanted, self.fake_equiv_4d_label_data)
 
     @skipif(no_nilearn)
-    def test_signal_extraction_4d_(self):
+    def test_signal_extraction_4d(self):
         self._test_4d_label([[-5.0652173913, -5.44565217391, 5.50543478261],
-                             [0, -2, .5],
-                             [-.3333333, -1, 2.5],
-                             [0, -2, .5],
-                             [-1.3333333, -5, 1]], self.fake_4d_label_data)
+                             [-7.02173913043, 11.1847826087, -4.33152173913],
+                             [-19.0869565217, 21.2391304348, -4.57608695652],
+                             [5.19565217391, -3.66304347826, -1.51630434783],
+                             [-12.0, 3., 0.5]], self.fake_4d_label_data)
 
     def _test_4d_label(self, wanted, fake_labels):
         # setup
@@ -71,7 +71,7 @@ def _test_4d_label(self, wanted, fake_labels):
 
         # run
         stats.SignalExtraction(in_file=self.filenames['in_file'],
-                               label_file=self.filenames['4d_label_file'],
+                               label_files=self.filenames['4d_label_file'],
                                class_labels=self.labels).run()
 
         self.assert_expected_output(wanted)

From 3883fe9b2e99838502e8e591299b707c024b0264 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Sun, 18 Sep 2016 22:22:42 +0000
Subject: [PATCH 17/29] make python3 happy

---
 nipype/algorithms/stats.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 0c2ee78d05..0ebe00904b 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -89,7 +89,7 @@ def _process_inputs(self):
                 masker = nl.NiftiMapsMasker(label_data)
                 n_labels = label_data.shape[3]
             else: # 3d file
-                if np.amax(label_data) > 1: # 3d label file
+                if np.amax(label_data.get_data()) > 1: # 3d label file
                     masker = nl.NiftiLabelsMasker(label_data)
                     # assuming consecutive positive integers for regions
                     n_labels = np.amax(label_data.get_data())

From 397ce99ac85aedde02d35cfb5bf1530c4ee269a7 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Sun, 18 Sep 2016 22:34:45 +0000
Subject: [PATCH 18/29] add global signal param

---
 nipype/algorithms/stats.py            | 3 +++
 nipype/algorithms/tests/test_stats.py | 2 +-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 0ebe00904b..941d2d0ce8 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -39,6 +39,9 @@ class SignalExtractionInputSpec(BaseInterfaceInputSpec):
     out_file = File('signals.tsv', usedefault=True, exists=False,
                     mandatory=False, desc='The name of the file to output to. '
                     'signals.tsv by default')
+    include_global = traits.Bool(False, usedefault=True, mandatory=False,
+                                 desc='If True, include an extra column '
+                                 'labeled "global"')
     detrend = traits.Bool(False, usedefault=True, mandatory=False,
                           desc='If True, perform detrending using nilearn.')
 
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index e198d087d5..f1282b53e5 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -36,7 +36,6 @@ def setUp(self):
         utils.save_toy_nii(self.fake_fmri_data, self.filenames['in_file'])
         utils.save_toy_nii(self.fake_label_data, self.filenames['label_files'])
 
-    @skipif(no_nilearn)
     def test_signal_extraction(self):
         # run
         stats.SignalExtraction(in_file=self.filenames['in_file'],
@@ -45,6 +44,7 @@ def test_signal_extraction(self):
         # assert
         self.assert_expected_output(self.base_wanted)
 
+
     @skipif(no_nilearn)
     @raises(ValueError)
     def test_signal_extraction_bad_label_list(self):

From 77f68cf595f2cbf6208ba66ccbe17d5ed5d43239 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Sun, 18 Sep 2016 22:40:02 +0000
Subject: [PATCH 19/29] add global signal test

---
 nipype/algorithms/tests/test_stats.py | 29 ++++++++++++++++++++++-----
 1 file changed, 24 insertions(+), 5 deletions(-)

diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index f1282b53e5..d8ed16bdd1 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -42,7 +42,7 @@ def test_signal_extraction(self):
                                label_files=self.filenames['label_files'],
                                class_labels=self.labels).run()
         # assert
-        self.assert_expected_output(self.base_wanted)
+        self.assert_expected_output(self.labels, self.base_wanted)
 
 
     @skipif(no_nilearn)
@@ -65,6 +65,25 @@ def test_signal_extraction_4d(self):
                              [5.19565217391, -3.66304347826, -1.51630434783],
                              [-12.0, 3., 0.5]], self.fake_4d_label_data)
 
+    @skipif(True)
+    @skipif(no_nilearn)
+    def test_signal_extraction_include_global(self):
+        # wanted
+        wanted_global = [[3./8], [-3./8], [1./8], [-7./8], [-9./8]]
+        for i, vals in enumerate(self.base_wanted):
+            wanted_global[i].extend(vals)
+        wanted_labels = ['global']
+        wanted_labels.extend(self.labels)
+
+        # run
+        stats.SignalExtraction(in_file=self.filenames['in_file'],
+                               label_files=self.filenames['label_files'],
+                               class_labels=self.labels,
+                               include_global=True).run()
+
+        # assert
+        self.assert_expected_output(wanted_labels, wanted_global)
+
     def _test_4d_label(self, wanted, fake_labels):
         # setup
         utils.save_toy_nii(fake_labels, self.filenames['4d_label_file'])
@@ -74,19 +93,19 @@ def _test_4d_label(self, wanted, fake_labels):
                                label_files=self.filenames['4d_label_file'],
                                class_labels=self.labels).run()
 
-        self.assert_expected_output(wanted)
+        self.assert_expected_output(self.labels, wanted)
 
-    def assert_expected_output(self, wanted):
+    def assert_expected_output(self, labels, wanted):
         with open(self.filenames['out_file'], 'r') as output:
             got = [line.split() for line in output]
             labels_got = got.pop(0) # remove header
-            assert_equal(labels_got, self.labels)
+            assert_equal(labels_got, labels)
             assert_equal(len(got), self.fake_fmri_data.shape[3],
                          'num rows and num volumes')
             # convert from string to float
             got = [[float(num) for num in row] for row in got]
             for i, time in enumerate(got):
-                assert_equal(len(self.labels), len(time))
+                assert_equal(len(labels), len(time))
                 for j, segment in enumerate(time):
                     assert_almost_equal(segment, wanted[i][j], decimal=1)
 

From 728609302bc0cfe751caf8c9c898adef5ba1cd55 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Sun, 18 Sep 2016 23:16:59 +0000
Subject: [PATCH 20/29] include global signal

---
 nipype/algorithms/stats.py            | 41 ++++++++++++++++++++-------
 nipype/algorithms/tests/test_stats.py |  1 -
 2 files changed, 30 insertions(+), 12 deletions(-)

diff --git a/nipype/algorithms/stats.py b/nipype/algorithms/stats.py
index 941d2d0ce8..7ff5c01040 100644
--- a/nipype/algorithms/stats.py
+++ b/nipype/algorithms/stats.py
@@ -62,15 +62,22 @@ class SignalExtraction(BaseInterface):
     >>> segments = ['CSF', 'gray', 'white']
     >>> seinterface.inputs.class_labels = segments
     >>> seinterface.inputs.detrend = True
+    >>> seinterface.inputs.include_global = True
     '''
     input_spec = SignalExtractionInputSpec
     output_spec = SignalExtractionOutputSpec
 
     def _run_interface(self, runtime):
-        masker = self._process_inputs()
+        masker, global_masker = self._process_inputs()
 
         region_signals = masker.fit_transform(self.inputs.in_file)
 
+        if global_masker:
+            self.inputs.class_labels.insert(0, 'global')
+            global_masker.fit()
+            global_signal= global_masker.transform(self.inputs.in_file)
+            region_signals = np.hstack((global_signal, region_signals))
+
         output = np.vstack((self.inputs.class_labels, region_signals.astype(str)))
 
         # save output
@@ -87,18 +94,19 @@ def _process_inputs(self):
             masker = nl.NiftiMapsMasker(self.inputs.label_files)
             n_labels = len(self.inputs.label_files)
         else: # list of size one, containing either a 3d or a 4d file
-            label_data = nb.load(self.inputs.label_files[0])
-            if len(label_data.shape) == 4: # 4d file
-                masker = nl.NiftiMapsMasker(label_data)
-                n_labels = label_data.shape[3]
+            self.label_data = nb.load(self.inputs.label_files[0])
+            if len(self.label_data.shape) == 4: # 4d file
+                masker = nl.NiftiMapsMasker(self.label_data)
+                n_labels = self.label_data.shape[3]
             else: # 3d file
-                if np.amax(label_data.get_data()) > 1: # 3d label file
-                    masker = nl.NiftiLabelsMasker(label_data)
+                if np.amax(self.label_data.get_data()) > 1: # 3d label file
+                    masker = nl.NiftiLabelsMasker(self.label_data)
                     # assuming consecutive positive integers for regions
-                    n_labels = np.amax(label_data.get_data())
+                    n_labels = np.amax(self.label_data.get_data())
                 else: # most probably a single probability map for one label
-                    masker = nl.NiftiMapsMasker(label_data)
+                    masker = nl.NiftiMapsMasker(self.label_data)
                     n_labels = 1
+        masker.set_params(detrend=self.inputs.detrend)
 
         # check label list size
         if len(self.inputs.class_labels) != n_labels:
@@ -108,8 +116,19 @@ def _process_inputs(self):
                                                     n_labels,
                                                     self.inputs.label_files))
 
-        masker.set_params(detrend=self.inputs.detrend)
-        return masker
+        if self.inputs.include_global:
+            all_ones_mask = nb.Nifti1Image(np.ones(self._label_data_shape()), np.eye(4))
+            global_masker = nl.NiftiLabelsMasker(all_ones_mask, detrend=self.inputs.detrend)
+        else:
+            global_masker = False
+
+        return masker, global_masker
+
+    def _label_data_shape(self):
+        if self.label_data:
+            return self.label_data.shape
+        else:
+            return nb.load(self.inputs.label_files[0]).shape
 
     def _list_outputs(self):
         outputs = self._outputs().get()
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index d8ed16bdd1..2d648f388d 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -65,7 +65,6 @@ def test_signal_extraction_4d(self):
                              [5.19565217391, -3.66304347826, -1.51630434783],
                              [-12.0, 3., 0.5]], self.fake_4d_label_data)
 
-    @skipif(True)
     @skipif(no_nilearn)
     def test_signal_extraction_include_global(self):
         # wanted

From 18c57cc2e444eff137faece55be7f5df40c752b1 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Mon, 19 Sep 2016 00:05:26 +0000
Subject: [PATCH 21/29] add back skipif no nilearn that disappeared

---
 nipype/algorithms/tests/test_stats.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/algorithms/tests/test_stats.py
index 2d648f388d..67e9438307 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/algorithms/tests/test_stats.py
@@ -36,6 +36,7 @@ def setUp(self):
         utils.save_toy_nii(self.fake_fmri_data, self.filenames['in_file'])
         utils.save_toy_nii(self.fake_label_data, self.filenames['label_files'])
 
+    @skipif(no_nilearn)
     def test_signal_extraction(self):
         # run
         stats.SignalExtraction(in_file=self.filenames['in_file'],

From 34fb1342d028fcc1884983d2a9b89945e6e0575e Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Mon, 19 Sep 2016 20:20:42 +0000
Subject: [PATCH 22/29] move to interfaces

---
 nipype/{algorithms/stats.py => interfaces/nilearn.py}  |  0
 .../test_stats.py => interfaces/tests/test_nilearn.py} | 10 +++++-----
 2 files changed, 5 insertions(+), 5 deletions(-)
 rename nipype/{algorithms/stats.py => interfaces/nilearn.py} (100%)
 rename nipype/{algorithms/tests/test_stats.py => interfaces/tests/test_nilearn.py} (95%)

diff --git a/nipype/algorithms/stats.py b/nipype/interfaces/nilearn.py
similarity index 100%
rename from nipype/algorithms/stats.py
rename to nipype/interfaces/nilearn.py
diff --git a/nipype/algorithms/tests/test_stats.py b/nipype/interfaces/tests/test_nilearn.py
similarity index 95%
rename from nipype/algorithms/tests/test_stats.py
rename to nipype/interfaces/tests/test_nilearn.py
index 67e9438307..bef4c60cf3 100644
--- a/nipype/algorithms/tests/test_stats.py
+++ b/nipype/interfaces/tests/test_nilearn.py
@@ -9,7 +9,7 @@
 
 from ...testing import (assert_equal, utils, assert_almost_equal, raises,
                         skipif)
-from .. import stats
+from .. import nilearn as iface
 
 no_nilearn = True
 try:
@@ -39,7 +39,7 @@ def setUp(self):
     @skipif(no_nilearn)
     def test_signal_extraction(self):
         # run
-        stats.SignalExtraction(in_file=self.filenames['in_file'],
+        iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['label_files'],
                                class_labels=self.labels).run()
         # assert
@@ -50,7 +50,7 @@ def test_signal_extraction(self):
     @raises(ValueError)
     def test_signal_extraction_bad_label_list(self):
         # run
-        stats.SignalExtraction(in_file=self.filenames['in_file'],
+        iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['label_files'],
                                class_labels=['bad']).run()
 
@@ -76,7 +76,7 @@ def test_signal_extraction_include_global(self):
         wanted_labels.extend(self.labels)
 
         # run
-        stats.SignalExtraction(in_file=self.filenames['in_file'],
+        iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['label_files'],
                                class_labels=self.labels,
                                include_global=True).run()
@@ -89,7 +89,7 @@ def _test_4d_label(self, wanted, fake_labels):
         utils.save_toy_nii(fake_labels, self.filenames['4d_label_file'])
 
         # run
-        stats.SignalExtraction(in_file=self.filenames['in_file'],
+        iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['4d_label_file'],
                                class_labels=self.labels).run()
 

From ab1e5f107f1ab4aa5b7c48255d56559e278c9ed3 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Tue, 20 Sep 2016 00:39:12 +0000
Subject: [PATCH 23/29] document include_global option more clearly

---
 nipype/interfaces/nilearn.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/nipype/interfaces/nilearn.py b/nipype/interfaces/nilearn.py
index 7ff5c01040..1b97573d93 100644
--- a/nipype/interfaces/nilearn.py
+++ b/nipype/interfaces/nilearn.py
@@ -41,7 +41,8 @@ class SignalExtractionInputSpec(BaseInterfaceInputSpec):
                     'signals.tsv by default')
     include_global = traits.Bool(False, usedefault=True, mandatory=False,
                                  desc='If True, include an extra column '
-                                 'labeled "global"')
+                                 'labeled "global", with values calculated from the entire brain '
+                                 '(instead of just regions)')
     detrend = traits.Bool(False, usedefault=True, mandatory=False,
                           desc='If True, perform detrending using nilearn.')
 

From 97dd3aace8274606fb0b434cc1bdc045a56b9acc Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Tue, 20 Sep 2016 06:19:27 +0000
Subject: [PATCH 24/29] actual global signal calculation

---
 nipype/interfaces/nilearn.py            | 44 +++++++++----------------
 nipype/interfaces/tests/test_nilearn.py | 40 ++++++++++++++--------
 2 files changed, 43 insertions(+), 41 deletions(-)

diff --git a/nipype/interfaces/nilearn.py b/nipype/interfaces/nilearn.py
index 1b97573d93..7fd8671c6c 100644
--- a/nipype/interfaces/nilearn.py
+++ b/nipype/interfaces/nilearn.py
@@ -73,7 +73,7 @@ def _run_interface(self, runtime):
 
         region_signals = masker.fit_transform(self.inputs.in_file)
 
-        if global_masker:
+        if global_masker != None:
             self.inputs.class_labels.insert(0, 'global')
             global_masker.fit()
             global_signal= global_masker.transform(self.inputs.in_file)
@@ -87,26 +87,19 @@ def _run_interface(self, runtime):
 
     def _process_inputs(self):
         ''' validate and  process inputs into useful form '''
-
         import nilearn.input_data as nl
+        import nilearn.image as nli
+
+        label_datas = [nb.load(nifti) for nifti in self.inputs.label_files]
+        label_data = nli.concat_imgs(label_datas)
 
         # determine form of label files, choose appropriate nilearn masker
-        if len(self.inputs.label_files) > 1: # list of 3D nifti images
-            masker = nl.NiftiMapsMasker(self.inputs.label_files)
-            n_labels = len(self.inputs.label_files)
-        else: # list of size one, containing either a 3d or a 4d file
-            self.label_data = nb.load(self.inputs.label_files[0])
-            if len(self.label_data.shape) == 4: # 4d file
-                masker = nl.NiftiMapsMasker(self.label_data)
-                n_labels = self.label_data.shape[3]
-            else: # 3d file
-                if np.amax(self.label_data.get_data()) > 1: # 3d label file
-                    masker = nl.NiftiLabelsMasker(self.label_data)
-                    # assuming consecutive positive integers for regions
-                    n_labels = np.amax(self.label_data.get_data())
-                else: # most probably a single probability map for one label
-                    masker = nl.NiftiMapsMasker(self.label_data)
-                    n_labels = 1
+        if len(label_datas) == 1 and np.amax(label_data.get_data()) > 1: # 3d label file
+            n_labels = np.amax(label_data.get_data())
+            masker = nl.NiftiLabelsMasker(label_data)
+        else: # one 4d file
+            n_labels = label_data.get_data().shape[3]
+            masker = nl.NiftiMapsMasker(label_data)
         masker.set_params(detrend=self.inputs.detrend)
 
         # check label list size
@@ -117,20 +110,15 @@ def _process_inputs(self):
                                                     n_labels,
                                                     self.inputs.label_files))
 
+        global_masker = None
         if self.inputs.include_global:
-            all_ones_mask = nb.Nifti1Image(np.ones(self._label_data_shape()), np.eye(4))
-            global_masker = nl.NiftiLabelsMasker(all_ones_mask, detrend=self.inputs.detrend)
-        else:
-            global_masker = False
+            global_label_data = label_data.get_data().clip(0, 1).sum(axis=3)
+            global_label_data = global_label_data[:, :, :, np.newaxis] # add back 4th dimension
+            global_label_data = nb.Nifti1Image(global_label_data, np.eye(4))
+            global_masker = nl.NiftiMapsMasker(global_label_data, detrend=self.inputs.detrend)
 
         return masker, global_masker
 
-    def _label_data_shape(self):
-        if self.label_data:
-            return self.label_data.shape
-        else:
-            return nb.load(self.inputs.label_files[0]).shape
-
     def _list_outputs(self):
         outputs = self._outputs().get()
         outputs['out_file'] = self.inputs.out_file
diff --git a/nipype/interfaces/tests/test_nilearn.py b/nipype/interfaces/tests/test_nilearn.py
index bef4c60cf3..6a3699c81e 100644
--- a/nipype/interfaces/tests/test_nilearn.py
+++ b/nipype/interfaces/tests/test_nilearn.py
@@ -27,6 +27,7 @@ class TestSignalExtraction(unittest.TestCase):
         'out_file': 'signals.tsv'
     }
     labels = ['csf', 'gray', 'white']
+    global_labels = ['global'] + labels
 
     def setUp(self):
         self.orig_dir = os.getcwd()
@@ -60,20 +61,14 @@ def test_signal_extraction_equiv_4d(self):
 
     @skipif(no_nilearn)
     def test_signal_extraction_4d(self):
-        self._test_4d_label([[-5.0652173913, -5.44565217391, 5.50543478261],
-                             [-7.02173913043, 11.1847826087, -4.33152173913],
-                             [-19.0869565217, 21.2391304348, -4.57608695652],
-                             [5.19565217391, -3.66304347826, -1.51630434783],
-                             [-12.0, 3., 0.5]], self.fake_4d_label_data)
+        self._test_4d_label(self.fourd_wanted, self.fake_4d_label_data)
 
     @skipif(no_nilearn)
-    def test_signal_extraction_include_global(self):
+    def test_signal_extraction_global(self):
         # wanted
-        wanted_global = [[3./8], [-3./8], [1./8], [-7./8], [-9./8]]
+        wanted_global = [[-4./6], [-1./6], [3./6], [-1./6], [-7./6]]
         for i, vals in enumerate(self.base_wanted):
             wanted_global[i].extend(vals)
-        wanted_labels = ['global']
-        wanted_labels.extend(self.labels)
 
         # run
         iface.SignalExtraction(in_file=self.filenames['in_file'],
@@ -82,18 +77,30 @@ def test_signal_extraction_include_global(self):
                                include_global=True).run()
 
         # assert
-        self.assert_expected_output(wanted_labels, wanted_global)
+        self.assert_expected_output(self.global_labels, wanted_global)
+
+    @skipif(no_nilearn)
+    def test_signal_extraction_4d_global(self):
+        # wanted
+        wanted_global = [[3./8], [-3./8], [1./8], [-7./8], [-9./8]]
+        for i, vals in enumerate(self.fourd_wanted):
+            wanted_global[i].extend(vals)
 
-    def _test_4d_label(self, wanted, fake_labels):
+        # run
+        self._test_4d_label(wanted_global, self.fake_4d_label_data, include_global=True)
+
+    def _test_4d_label(self, wanted, fake_labels, include_global=False):
         # setup
         utils.save_toy_nii(fake_labels, self.filenames['4d_label_file'])
 
         # run
         iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['4d_label_file'],
-                               class_labels=self.labels).run()
+                               class_labels=self.labels,
+                               include_global=include_global).run()
 
-        self.assert_expected_output(self.labels, wanted)
+        wanted_labels = self.global_labels if include_global else self.labels
+        self.assert_expected_output(wanted_labels, wanted)
 
     def assert_expected_output(self, labels, wanted):
         with open(self.filenames['out_file'], 'r') as output:
@@ -159,3 +166,10 @@ def tearDown(self):
 
                                     [[0.3, 0.3, 0.4],
                                      [0.3, 0.4, 0.3]]]])
+
+
+    fourd_wanted = [[-5.0652173913, -5.44565217391, 5.50543478261],
+                    [-7.02173913043, 11.1847826087, -4.33152173913],
+                    [-19.0869565217, 21.2391304348, -4.57608695652],
+                    [5.19565217391, -3.66304347826, -1.51630434783],
+                    [-12.0, 3., 0.5]]

From 74b804e290bfab21b44f87a39edb280f97bfdc5b Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Tue, 20 Sep 2016 19:41:47 +0000
Subject: [PATCH 25/29] add new param incl_shared_variance

---
 nipype/interfaces/nilearn.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/nipype/interfaces/nilearn.py b/nipype/interfaces/nilearn.py
index 7fd8671c6c..72d27c948c 100644
--- a/nipype/interfaces/nilearn.py
+++ b/nipype/interfaces/nilearn.py
@@ -39,6 +39,11 @@ class SignalExtractionInputSpec(BaseInterfaceInputSpec):
     out_file = File('signals.tsv', usedefault=True, exists=False,
                     mandatory=False, desc='The name of the file to output to. '
                     'signals.tsv by default')
+    incl_shared_variance = traits.Bool(True, usedefault=True, mandatory=False, desc='By default '
+                                       '(True), returns simple time series calculated from each '
+                                       'region independently (e.g., for noise regression). If '
+                                       'False, returns unique signals for each region, discarding '
+                                       'shared variance (e.g., for connectivity)')
     include_global = traits.Bool(False, usedefault=True, mandatory=False,
                                  desc='If True, include an extra column '
                                  'labeled "global", with values calculated from the entire brain '

From 5acb6dc4f4951ebcce016c2061c1c8ff4cb036d4 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Tue, 20 Sep 2016 20:25:35 +0000
Subject: [PATCH 26/29] modify tests for new param

---
 nipype/interfaces/nilearn.py            |  3 +-
 nipype/interfaces/tests/test_nilearn.py | 58 +++++++++++++++++--------
 2 files changed, 43 insertions(+), 18 deletions(-)

diff --git a/nipype/interfaces/nilearn.py b/nipype/interfaces/nilearn.py
index 72d27c948c..a2e1593169 100644
--- a/nipype/interfaces/nilearn.py
+++ b/nipype/interfaces/nilearn.py
@@ -47,7 +47,8 @@ class SignalExtractionInputSpec(BaseInterfaceInputSpec):
     include_global = traits.Bool(False, usedefault=True, mandatory=False,
                                  desc='If True, include an extra column '
                                  'labeled "global", with values calculated from the entire brain '
-                                 '(instead of just regions)')
+                                 '(instead of just regions). Only has effect with 4D probability '
+                                 'maps.')
     detrend = traits.Bool(False, usedefault=True, mandatory=False,
                           desc='If True, perform detrending using nilearn.')
 
diff --git a/nipype/interfaces/tests/test_nilearn.py b/nipype/interfaces/tests/test_nilearn.py
index 6a3699c81e..8bf200a58c 100644
--- a/nipype/interfaces/tests/test_nilearn.py
+++ b/nipype/interfaces/tests/test_nilearn.py
@@ -38,34 +38,38 @@ def setUp(self):
         utils.save_toy_nii(self.fake_label_data, self.filenames['label_files'])
 
     @skipif(no_nilearn)
-    def test_signal_extraction(self):
+    def test_signal_extract_no_shared(self):
         # run
         iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['label_files'],
-                               class_labels=self.labels).run()
+                               class_labels=self.labels,
+                               incl_shared_variance=False).run()
         # assert
         self.assert_expected_output(self.labels, self.base_wanted)
 
 
     @skipif(no_nilearn)
     @raises(ValueError)
-    def test_signal_extraction_bad_label_list(self):
+    def test_signal_extr_bad_label_list(self):
         # run
         iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['label_files'],
-                               class_labels=['bad']).run()
+                               class_labels=['bad'],
+                               incl_shared_variance=False).run()
 
     @skipif(no_nilearn)
-    def test_signal_extraction_equiv_4d(self):
-        self._test_4d_label(self.base_wanted, self.fake_equiv_4d_label_data)
+    def test_signal_extr_equiv_4d_no_shared(self):
+        self._test_4d_label(self.base_wanted, self.fake_equiv_4d_label_data,
+                            incl_shared_variance=False)
 
     @skipif(no_nilearn)
-    def test_signal_extraction_4d(self):
-        self._test_4d_label(self.fourd_wanted, self.fake_4d_label_data)
+    def test_signal_extr_4d_no_shared(self):
+        # set up & run & assert
+        self._test_4d_label(self.fourd_wanted, self.fake_4d_label_data, incl_shared_variance=False)
 
     @skipif(no_nilearn)
-    def test_signal_extraction_global(self):
-        # wanted
+    def test_signal_extr_global_no_shared(self):
+        # set up
         wanted_global = [[-4./6], [-1./6], [3./6], [-1./6], [-7./6]]
         for i, vals in enumerate(self.base_wanted):
             wanted_global[i].extend(vals)
@@ -74,32 +78,52 @@ def test_signal_extraction_global(self):
         iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['label_files'],
                                class_labels=self.labels,
-                               include_global=True).run()
+                               include_global=True,
+                               incl_shared_variance=False).run()
 
         # assert
         self.assert_expected_output(self.global_labels, wanted_global)
 
     @skipif(no_nilearn)
-    def test_signal_extraction_4d_global(self):
-        # wanted
+    def test_signal_extr_4d_global_no_shared(self):
+        # set up
         wanted_global = [[3./8], [-3./8], [1./8], [-7./8], [-9./8]]
         for i, vals in enumerate(self.fourd_wanted):
             wanted_global[i].extend(vals)
 
-        # run
-        self._test_4d_label(wanted_global, self.fake_4d_label_data, include_global=True)
+        # run & assert
+        self._test_4d_label(wanted_global, self.fake_4d_label_data,
+                            include_global=True, incl_shared_variance=False)
 
-    def _test_4d_label(self, wanted, fake_labels, include_global=False):
-        # setup
+    @skipif(no_nilearn)
+    def test_signal_extr_shared(self):
+        # set up
+        wanted = []
+        for vol in range(self.fake_fmri_data.shape[3]):
+            volume = self.fake_fmri_data[:, :, :, vol].flatten()
+            wanted_row = []
+            for reg in range(self.fake_4d_label_data.shape[3]):
+                region = self.fake_4d_label_data[:, :, :, reg].flatten()
+                wanted_row.append(np.average(volume, weights=region))
+            wanted.append(wanted_row)
+
+        # run & assert
+        self._test_4d_label(wanted, self.fake_4d_label_data)
+
+    def _test_4d_label(self, wanted, fake_labels, include_global=False, incl_shared_variance=True):
+        # set up
         utils.save_toy_nii(fake_labels, self.filenames['4d_label_file'])
 
         # run
         iface.SignalExtraction(in_file=self.filenames['in_file'],
                                label_files=self.filenames['4d_label_file'],
                                class_labels=self.labels,
+                               incl_shared_variance=incl_shared_variance,
                                include_global=include_global).run()
 
         wanted_labels = self.global_labels if include_global else self.labels
+
+        # assert
         self.assert_expected_output(wanted_labels, wanted)
 
     def assert_expected_output(self, labels, wanted):

From 9df8d8ae1df772723aa7e903fa38341230f01191 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Wed, 21 Sep 2016 04:19:43 +0000
Subject: [PATCH 27/29] calculate linear regression individually be default

---
 nipype/interfaces/nilearn.py            | 42 ++++++++++++++-----------
 nipype/interfaces/tests/test_nilearn.py |  4 +--
 2 files changed, 26 insertions(+), 20 deletions(-)

diff --git a/nipype/interfaces/nilearn.py b/nipype/interfaces/nilearn.py
index a2e1593169..0b54e95952 100644
--- a/nipype/interfaces/nilearn.py
+++ b/nipype/interfaces/nilearn.py
@@ -75,15 +75,12 @@ class SignalExtraction(BaseInterface):
     output_spec = SignalExtractionOutputSpec
 
     def _run_interface(self, runtime):
-        masker, global_masker = self._process_inputs()
+        maskers = self._process_inputs()
 
-        region_signals = masker.fit_transform(self.inputs.in_file)
-
-        if global_masker != None:
-            self.inputs.class_labels.insert(0, 'global')
-            global_masker.fit()
-            global_signal= global_masker.transform(self.inputs.in_file)
-            region_signals = np.hstack((global_signal, region_signals))
+        signals = []
+        for masker in maskers:
+            signals.append(masker.fit_transform(self.inputs.in_file))
+        region_signals = np.hstack(signals)
 
         output = np.vstack((self.inputs.class_labels, region_signals.astype(str)))
 
@@ -92,21 +89,26 @@ def _run_interface(self, runtime):
         return runtime
 
     def _process_inputs(self):
-        ''' validate and  process inputs into useful form '''
+        ''' validate and  process inputs into useful form.
+        Returns a list of nilearn maskers and the list of corresponding label names.'''
         import nilearn.input_data as nl
         import nilearn.image as nli
 
-        label_datas = [nb.load(nifti) for nifti in self.inputs.label_files]
-        label_data = nli.concat_imgs(label_datas)
+        label_data = nli.concat_imgs(self.inputs.label_files)
+        maskers = []
 
         # determine form of label files, choose appropriate nilearn masker
-        if len(label_datas) == 1 and np.amax(label_data.get_data()) > 1: # 3d label file
+        if np.amax(label_data.get_data()) > 1: # 3d label file
             n_labels = np.amax(label_data.get_data())
-            masker = nl.NiftiLabelsMasker(label_data)
-        else: # one 4d file
+            maskers.append(nl.NiftiLabelsMasker(label_data))
+        else: # 4d labels
             n_labels = label_data.get_data().shape[3]
-            masker = nl.NiftiMapsMasker(label_data)
-        masker.set_params(detrend=self.inputs.detrend)
+            if self.inputs.incl_shared_variance: # 4d labels, independent computation
+                for img in nli.iter_img(label_data):
+                    sortof_4d_img = nb.Nifti1Image(img.get_data()[:, :, :, np.newaxis], np.eye(4))
+                    maskers.append(nl.NiftiMapsMasker(sortof_4d_img))
+            else: # 4d labels, one computation fitting all
+                maskers.append(nl.NiftiMapsMasker(label_data))
 
         # check label list size
         if len(self.inputs.class_labels) != n_labels:
@@ -116,14 +118,18 @@ def _process_inputs(self):
                                                     n_labels,
                                                     self.inputs.label_files))
 
-        global_masker = None
         if self.inputs.include_global:
             global_label_data = label_data.get_data().clip(0, 1).sum(axis=3)
             global_label_data = global_label_data[:, :, :, np.newaxis] # add back 4th dimension
             global_label_data = nb.Nifti1Image(global_label_data, np.eye(4))
             global_masker = nl.NiftiMapsMasker(global_label_data, detrend=self.inputs.detrend)
+            maskers.insert(0, global_masker)
+            self.inputs.class_labels.insert(0, 'global')
+
+        for masker in maskers:
+            masker.set_params(detrend=self.inputs.detrend)
 
-        return masker, global_masker
+        return maskers
 
     def _list_outputs(self):
         outputs = self._outputs().get()
diff --git a/nipype/interfaces/tests/test_nilearn.py b/nipype/interfaces/tests/test_nilearn.py
index 8bf200a58c..074e71fee9 100644
--- a/nipype/interfaces/tests/test_nilearn.py
+++ b/nipype/interfaces/tests/test_nilearn.py
@@ -104,9 +104,9 @@ def test_signal_extr_shared(self):
             wanted_row = []
             for reg in range(self.fake_4d_label_data.shape[3]):
                 region = self.fake_4d_label_data[:, :, :, reg].flatten()
-                wanted_row.append(np.average(volume, weights=region))
-            wanted.append(wanted_row)
+                wanted_row.append((volume*region).sum()/(region*region).sum())
 
+            wanted.append(wanted_row)
         # run & assert
         self._test_4d_label(wanted, self.fake_4d_label_data)
 

From 9dad02e01bd800c5fa6d69529dff79027624552b Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Wed, 21 Sep 2016 21:10:52 +0000
Subject: [PATCH 28/29] binarize global maps

---
 nipype/interfaces/nilearn.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/nipype/interfaces/nilearn.py b/nipype/interfaces/nilearn.py
index 0b54e95952..19bd890245 100644
--- a/nipype/interfaces/nilearn.py
+++ b/nipype/interfaces/nilearn.py
@@ -119,7 +119,8 @@ def _process_inputs(self):
                                                     self.inputs.label_files))
 
         if self.inputs.include_global:
-            global_label_data = label_data.get_data().clip(0, 1).sum(axis=3)
+            global_label_data = label_data.get_data().sum(axis=3)
+            global_label_data = np.rint(global_label_data).astype(int).clip(0, 1)
             global_label_data = global_label_data[:, :, :, np.newaxis] # add back 4th dimension
             global_label_data = nb.Nifti1Image(global_label_data, np.eye(4))
             global_masker = nl.NiftiMapsMasker(global_label_data, detrend=self.inputs.detrend)

From 850ed4173d63cf6bd8688502fa8e8859fbdcc246 Mon Sep 17 00:00:00 2001
From: Shoshana Berleant <berleant@stanford.edu>
Date: Thu, 22 Sep 2016 01:24:40 +0000
Subject: [PATCH 29/29] clearer 4Dification, doc edit

---
 nipype/interfaces/nilearn.py | 35 +++++++++++++++++++----------------
 1 file changed, 19 insertions(+), 16 deletions(-)

diff --git a/nipype/interfaces/nilearn.py b/nipype/interfaces/nilearn.py
index 19bd890245..bc662826d0 100644
--- a/nipype/interfaces/nilearn.py
+++ b/nipype/interfaces/nilearn.py
@@ -25,10 +25,10 @@
 class SignalExtractionInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='4-D fMRI nii file')
     label_files = InputMultiPath(File(exists=True), mandatory=True,
-                                desc='a 3-D label image, with 0 denoting '
-                                'background, or a list of 3-D probability '
-                                'maps (one per label) or the equivalent 4D '
-                                'file.')
+                                 desc='a 3-D label image, with 0 denoting '
+                                 'background, or a list of 3-D probability '
+                                 'maps (one per label) or the equivalent 4D '
+                                 'file.')
     class_labels = traits.List(mandatory=True,
                                desc='Human-readable labels for each segment '
                                'in the label file, in order. The length of '
@@ -43,12 +43,12 @@ class SignalExtractionInputSpec(BaseInterfaceInputSpec):
                                        '(True), returns simple time series calculated from each '
                                        'region independently (e.g., for noise regression). If '
                                        'False, returns unique signals for each region, discarding '
-                                       'shared variance (e.g., for connectivity)')
+                                       'shared variance (e.g., for connectivity. Only has effect '
+                                       'with 4D probability maps.')
     include_global = traits.Bool(False, usedefault=True, mandatory=False,
                                  desc='If True, include an extra column '
                                  'labeled "global", with values calculated from the entire brain '
-                                 '(instead of just regions). Only has effect with 4D probability '
-                                 'maps.')
+                                 '(instead of just regions).')
     detrend = traits.Bool(False, usedefault=True, mandatory=False,
                           desc='If True, perform detrending using nilearn.')
 
@@ -105,8 +105,7 @@ def _process_inputs(self):
             n_labels = label_data.get_data().shape[3]
             if self.inputs.incl_shared_variance: # 4d labels, independent computation
                 for img in nli.iter_img(label_data):
-                    sortof_4d_img = nb.Nifti1Image(img.get_data()[:, :, :, np.newaxis], np.eye(4))
-                    maskers.append(nl.NiftiMapsMasker(sortof_4d_img))
+                    maskers.append(nl.NiftiMapsMasker(self._4d(img.get_data(), img.affine)))
             else: # 4d labels, one computation fitting all
                 maskers.append(nl.NiftiMapsMasker(label_data))
 
@@ -115,15 +114,14 @@ def _process_inputs(self):
             raise ValueError('The length of class_labels {} does not '
                              'match the number of regions {} found in '
                              'label_files {}'.format(self.inputs.class_labels,
-                                                    n_labels,
-                                                    self.inputs.label_files))
+                                                     n_labels,
+                                                     self.inputs.label_files))
 
         if self.inputs.include_global:
-            global_label_data = label_data.get_data().sum(axis=3)
-            global_label_data = np.rint(global_label_data).astype(int).clip(0, 1)
-            global_label_data = global_label_data[:, :, :, np.newaxis] # add back 4th dimension
-            global_label_data = nb.Nifti1Image(global_label_data, np.eye(4))
-            global_masker = nl.NiftiMapsMasker(global_label_data, detrend=self.inputs.detrend)
+            global_label_data = label_data.get_data().sum(axis=3) # sum across all regions
+            global_label_data = np.rint(global_label_data).astype(int).clip(0, 1) # binarize
+            global_label_data = self._4d(global_label_data, label_data.affine)
+            global_masker = nl.NiftiLabelsMasker(global_label_data, detrend=self.inputs.detrend)
             maskers.insert(0, global_masker)
             self.inputs.class_labels.insert(0, 'global')
 
@@ -132,6 +130,11 @@ def _process_inputs(self):
 
         return maskers
 
+    def _4d(self, array, affine):
+        ''' takes a 3-dimensional numpy array and an affine,
+        returns the equivalent 4th dimensional nifti file '''
+        return nb.Nifti1Image(array[:, :, :, np.newaxis], affine)
+
     def _list_outputs(self):
         outputs = self._outputs().get()
         outputs['out_file'] = self.inputs.out_file