Repository: bids-standard/pybids

Issue:
Convolve can return multiple regressors (not just one)
currently the implementation of Convolve assumes the hrf function returns an array with only one axis expected to have values. but if 'glover + derivatives' is included, I would expect two columns to be created (orig and orig_derivative). 

Discussion:

Patch:
diff --git a/bids/modeling/transformations/base.py b/bids/modeling/transformations/base.py
index 55588098..3c87f69f 100644
--- a/bids/modeling/transformations/base.py
+++ b/bids/modeling/transformations/base.py
@@ -313,8 +313,17 @@ class Transformation(metaclass=ABCMeta):
                 if self.output_suffix is not None:
                     _output += self.output_suffix
 
-                col.name = _output
-                self.collection[_output] = col
+                # If multiple variables were returned, add each one separately
+                if isinstance(result, (list, tuple)):
+                    # rename first output
+                    result[0].name = _output
+                    self.collection[_output] = result[0]
+
+                    for r in result[1:]:
+                        self.collection[r.name] = r
+                else:
+                    col.name = _output
+                    self.collection[_output] = col
 
     @abstractmethod
     def _transform(self, **kwargs):
diff --git a/bids/modeling/transformations/compute.py b/bids/modeling/transformations/compute.py
index 9478f926..71877c06 100644
--- a/bids/modeling/transformations/compute.py
+++ b/bids/modeling/transformations/compute.py
@@ -39,7 +39,7 @@ class Convolve(Transformation):
     -----
     Uses the HRF convolution functions implemented in nistats.
     """
-
+    _groupable = False
     _input_type = 'variable'
     _return_type = 'variable'
 
@@ -93,9 +93,16 @@ class Convolve(Transformation):
             oversampling=np.ceil(effective_sr / sampling_rate)
             )
 
-        return DenseRunVariable(
-            name=var.name, values=convolved[0], run_info=var.run_info,
-            source=var.source, sampling_rate=sampling_rate)
+        results = []
+        arr, names = convolved
+        for conv, name in zip(np.split(arr, arr.shape[1], axis=1), names):
+            new_name = '_'.join([var.name, name.split('_')[-1]]) if '_' in name else var.name
+            results.append(
+                DenseRunVariable(
+                    name=new_name, values=conv, run_info=var.run_info,
+                    source=var.source, sampling_rate=sampling_rate)
+            )
+        return results
 
 
 class Demean(Transformation):


Test patch:
diff --git a/bids/modeling/tests/test_transformations.py b/bids/modeling/tests/test_transformations.py
index e79ee6a0..91b15506 100644
--- a/bids/modeling/tests/test_transformations.py
+++ b/bids/modeling/tests/test_transformations.py
@@ -58,18 +58,22 @@ def sparse_run_variable_with_missing_values():
 
 def test_convolve(collection):
     rt = collection.variables['RT']
-    transform.Convolve(collection, 'RT', output='reaction_time')
+    transform.Convolve(collection, ['RT'], output=['reaction_time'])
     rt_conv = collection.variables['reaction_time']
 
     assert rt_conv.values.shape[0] == \
         rt.get_duration() * collection.sampling_rate
 
-    transform.ToDense(collection, 'RT', output='rt_dense')
-    transform.Convolve(collection, 'rt_dense', output='dense_convolved')
+    transform.ToDense(collection, ['RT'], output=['rt_dense'])
+    transform.Convolve(collection, 'rt_dense', derivative=True)
+
+    # test the derivative exists
+    assert collection.variables.get('rt_dense_derivative')
 
     dense_conv = collection.variables['reaction_time']
 
     assert dense_conv.values.shape[0] == \
+        collection.variables['rt_dense_derivative'].values.shape[0] == \
         rt.get_duration() * collection.sampling_rate
 
     # Test adapative oversampling computation
@@ -77,34 +81,40 @@ def test_convolve(collection):
     # To resolve 1Hz frequencies, we must sample at >=2Hz
     args = (mock.ANY, 'spm', mock.ANY)
     kwargs = dict(fir_delays=None, min_onset=0)
-    with mock.patch('bids.modeling.transformations.compute.hrf') as mocked:
+    mock_return = (np.array([[0, 0]]), ["cond"])
+    with mock.patch('bids.modeling.transformations.compute.hrf.compute_regressor') as mocked:
+        mocked.return_value = mock_return
         # Sampling rate is 10Hz, no oversampling needed
         transform.Convolve(collection, 'RT', output='rt_mock')
-        mocked.compute_regressor.assert_called_with(*args, oversampling=1.0, **kwargs)
+        mocked.assert_called_with(*args, oversampling=1.0, **kwargs)
 
-    with mock.patch('bids.modeling.transformations.compute.hrf') as mocked:
+    with mock.patch('bids.modeling.transformations.compute.hrf.compute_regressor') as mocked:
+        mocked.return_value = mock_return
         # Sampling rate is 10Hz, no oversampling needed
         transform.Convolve(collection, 'rt_dense', output='rt_mock')
-        mocked.compute_regressor.assert_called_with(*args, oversampling=1.0, **kwargs)
+        mocked.assert_called_with(*args, oversampling=1.0, **kwargs)
 
-    with mock.patch('bids.modeling.transformations.compute.hrf') as mocked:
+    with mock.patch('bids.modeling.transformations.compute.hrf.compute_regressor') as mocked:
+        mocked.return_value = mock_return
         # Slow sampling rate, oversample (4x) to 2Hz
         collection.sampling_rate = 0.5
         transform.Convolve(collection, 'RT', output='rt_mock')
-        mocked.compute_regressor.assert_called_with(*args, oversampling=4.0, **kwargs)
+        mocked.assert_called_with(*args, oversampling=4.0, **kwargs)
 
-    with mock.patch('bids.modeling.transformations.compute.hrf') as mocked:
+    with mock.patch('bids.modeling.transformations.compute.hrf.compute_regressor') as mocked:
+        mocked.return_value = mock_return
         # Dense variable is already sampled at 10Hz, no oversampling needed
         collection.sampling_rate = 0.5
         transform.Convolve(collection, 'rt_dense', output='rt_mock')
-        mocked.compute_regressor.assert_called_with(*args, oversampling=1.0, **kwargs)
+        mocked.assert_called_with(*args, oversampling=1.0, **kwargs)
 
-    with mock.patch('bids.modeling.transformations.compute.hrf') as mocked:
+    with mock.patch('bids.modeling.transformations.compute.hrf.compute_regressor') as mocked:
+        mocked.return_value = mock_return
         # Onset requires 10Hz resolution, oversample (2x) to 20Hz
         collection.sampling_rate = 10
         collection['RT'].onset[0] += 0.1
         transform.Convolve(collection, 'RT', output='rt_mock')
-        mocked.compute_regressor.assert_called_with(*args, oversampling=2.0, **kwargs)
+        mocked.assert_called_with(*args, oversampling=2.0, **kwargs)
 
 
 def test_convolve_impulse():
@@ -126,7 +136,7 @@ def test_convolve_impulse():
 def test_rename(collection):
     dense_rt = collection.variables['RT'].to_dense(collection.sampling_rate)
     assert len(dense_rt.values) == math.ceil(len(SUBJECTS) * NRUNS * SCAN_LENGTH * collection.sampling_rate)
-    transform.Rename(collection, 'RT', output='reaction_time')
+    transform.Rename(collection, ['RT'], output=['reaction_time'])
     assert 'reaction_time' in collection.variables
     assert 'RT' not in collection.variables
     col = collection.variables['reaction_time']