ENH: Enable models for sparsely sampled fMRI series

bids/analysis/analysis.py

@@ -388,18 +388,29 @@ def get_design_matrix(self, names=None, format='long', mode='both',
 
             if sampling_rate == 'TR':
                 trs = {var.run_info[0].tr for var in self.collection.variables.values()}
+                tas = {var.run_info[0].ta for var in self.collection.variables.values()}
                 if not trs:
                     raise ValueError("Repetition time unavailable; specify sampling_rate "
                                      "explicitly")
                 elif len(trs) > 1:
                     raise ValueError("Non-unique Repetition times found ({!r}); specify "
-                                     "sampling_rate explicitly")
-                sampling_rate = 1. / trs.pop()
+                                     "sampling_rate explicitly".format(trs))
+                TR = trs.pop()
+                if not tas:
+                    warnings.warn("Acquisition time unavailable; assuming TA = TR")
+                    tas = {TR}
+                elif len(tas) > 1:
+                    raise ValueError("Non-unique acquisition times found ({!r})".format(tas))
+
+                sampling_rate = 1. / TR
+                acquisition_time = tas.pop()
             elif sampling_rate == 'highest':
                 sampling_rate = None
             dense_df = coll.to_df(names, format='wide',
                                   include_sparse=include_sparse,
-                                  sampling_rate=sampling_rate, **kwargs)
+                                  sampling_rate=sampling_rate,
+                                  integration_window=acquisition_time,
+                                  **kwargs)
             if dense_df is not None:
                 dense_df = dense_df.drop(['onset', 'duration'], axis=1)
 

bids/analysis/tests/test_transformations.py

@@ -28,7 +28,7 @@ def sparse_run_variable_with_missing_values():
         'duration': [1.2, 1.6, 0.8, 2],
         'amplitude': [1, 1, np.nan, 1]
     })
-    run_info = [RunInfo({'subject': '01'}, 20, 2, 'dummy.nii.gz')]
+    run_info = [RunInfo({'subject': '01'}, 20, 2, 2, 'dummy.nii.gz')]
     var = SparseRunVariable('var', data, run_info, 'events')
     return BIDSRunVariableCollection([var])
 

bids/analysis/transformations/compute.py

@@ -2,6 +2,7 @@
 Transformations that primarily involve numerical computation on variables.
 '''
 
+from math import gcd
 import numpy as np
 import pandas as pd
 from bids.utils import listify
@@ -35,6 +36,25 @@ def _transform(self, var, model='spm', derivative=False, dispersion=False,
 
         if isinstance(var, SparseRunVariable):
             sr = self.collection.sampling_rate
+            # Resolve diferences in TR and TA to milliseconds
+            trs = {ri.tr for ri in var.run_info}
+            tas = {ri.ta for ri in var.run_info}
+            assert len(trs) == 1
+            assert len(tas) == 1
+            TR = int(np.round(1000. * trs.pop()))
+            TA = int(np.round(1000. * tas.pop()))
+            SR = int(np.round(1000. / sr))
+            if TA is None or TA < TR:
+                # Use a unit that fits an whole number of times into both
+                # the interscan interval (TR) and the integration window (TA)
+                dt = gcd(TR, TA)
+                if dt > SR:
+                    # Find the nearest whole factor of dt >= SR
+                    # This compromises on the target sampling rate and
+                    # one that neatly bins TR and TA
+                    dt = dt // (dt // SR)
+                sr = 1000. / dt
+
             var = var.to_dense(sr)
 
         df = var.to_df(entities=False)

bids/variables/entities.py

@@ -36,23 +36,25 @@ class RunNode(Node):
         image_file (str): The full path to the corresponding nifti image.
         duration (float): Duration of the run, in seconds.
         repetition_time (float): TR for the run.
+        acquisition_time (float): TA for the run.
         task (str): The task name for this run.
     '''
 
-    def __init__(self, entities, image_file, duration, repetition_time):
+    def __init__(self, entities, image_file, duration, repetition_time, acquisition_time=None):
         self.image_file = image_file
         self.duration = duration
         self.repetition_time = repetition_time
+        self.acquisition_time = acquisition_time or repetition_time
         super(RunNode, self).__init__('run', entities)
 
     def get_info(self):
 
         return RunInfo(self.entities, self.duration, self.repetition_time,
-                       self.image_file)
+                       self.acquisition_time, self.image_file)
 
 
 # Stores key information for each Run.
-RunInfo = namedtuple('RunInfo', ['entities', 'duration', 'tr', 'image'])
+RunInfo = namedtuple('RunInfo', ['entities', 'duration', 'tr', 'ta', 'image'])
 
 
 class NodeIndex(Node):

bids/variables/io.py

@@ -82,6 +82,27 @@ def load_variables(layout, types=None, levels=None, skip_empty=True,
     return dataset
 
 
+def _get_timing_info(img_md):
+    if 'RepetitionTime' in img_md:
+        tr = img_md['RepetitionTime']
+        if 'DelayTime' in img_md:
+            ta = tr - img_md['DelayTime']
+        elif 'SliceTiming' in img_md:
+            slicetimes = sorted(img_md['SliceTiming'])
+            # a, b ... z
+            # z = final slice onset, b - a = slice duration
+            ta = np.round(slicetimes[-1] + slicetimes[1] - slicetimes[0], 3)
+            # If the "silence" is <1ms, consider acquisition continuous
+            if np.abs(tr - ta) < 1e-3:
+                ta = tr
+        else:
+            ta = tr
+    elif 'VolumeTiming' in img_md:
+        return NotImplemented
+
+    return tr, ta
+
+
 def _load_time_variables(layout, dataset=None, columns=None, scan_length=None,
                          drop_na=True, events=True, physio=True, stim=True,
                          regressors=True, skip_empty=True, **selectors):
@@ -141,8 +162,9 @@ def _load_time_variables(layout, dataset=None, columns=None, scan_length=None,
         if 'run' in entities:
             entities['run'] = int(entities['run'])
 
-        tr = layout.get_metadata(img_f, suffix='bold', domains=domains,
-                                 full_search=True)['RepetitionTime']
+        img_md = layout.get_metadata(img_f, suffix='bold', domains=domains,
+                                     full_search=True)
+        tr, ta = _get_timing_info(img_md)
 
         # Get duration of run: first try to get it directly from the image
         # header; if that fails, try to get NumberOfVolumes from the
@@ -162,7 +184,8 @@ def _load_time_variables(layout, dataset=None, columns=None, scan_length=None,
                 raise ValueError(msg)
 
         run = dataset.get_or_create_node('run', entities, image_file=img_f,
-                                         duration=duration, repetition_time=tr)
+                                         duration=duration, repetition_time=tr,
+                                         acquisition_time=ta)
         run_info = run.get_info()
 
         # Process event files

bids/variables/kollekshuns.py

@@ -243,7 +243,7 @@ def _all_dense(self):
                     for v in self.variables.values()])
 
     def resample(self, sampling_rate=None, variables=None, force_dense=False,
-                 in_place=False, kind='linear'):
+                 in_place=False, kind='linear', **kwargs):
         ''' Resample all dense variables (and optionally, sparse ones) to the
         specified sampling rate.
 
@@ -276,7 +276,8 @@ def resample(self, sampling_rate=None, variables=None, force_dense=False,
                 # None if in_place; no update needed
                 _var = var.resample(sampling_rate,
                                     inplace=in_place,
-                                    kind=kind)
+                                    kind=kind,
+                                    **kwargs)
                 if not in_place:
                     _variables[name] = _var
 
@@ -289,6 +290,7 @@ def resample(self, sampling_rate=None, variables=None, force_dense=False,
 
     def to_df(self, variables=None, format='wide', sparse=True,
               sampling_rate=None, include_sparse=True, include_dense=True,
+              integration_window=None,
               **kwargs):
         ''' Merge columns into a single pandas DataFrame.
 
@@ -343,9 +345,10 @@ def to_df(self, variables=None, format='wide', sparse=True,
             sampling_rate = sampling_rate or self.sampling_rate
 
             # Make sure all variables have the same sampling rate
-            variables = list(self.resample(sampling_rate, variables,
-                                           force_dense=True,
-                                           in_place=False).values())
+            variables = list(
+                self.resample(sampling_rate, variables,
+                              force_dense=True, in_place=False,
+                              integration_window=integration_window).values())
 
         return super(BIDSRunVariableCollection, self).to_df(variables, format,
                                                             **kwargs)

bids/variables/tests/test_variables.py

@@ -19,7 +19,7 @@ def generate_DEV(name='test', sr=20, duration=480):
     ent_names = ['task', 'run', 'session', 'subject']
     entities = {e: uuid.uuid4().hex for e in ent_names}
     image = uuid.uuid4().hex + '.nii.gz'
-    run_info = RunInfo(entities, duration, 2, image)
+    run_info = RunInfo(entities, duration, 2, 2, image)
     return DenseRunVariable('test', values, run_info, 'dummy', sr)
 
 

bids/variables/variables.py

@@ -253,8 +253,8 @@ def split(self, grouper):
         subsets = []
         for i, (name, g) in enumerate(data.groupby(grouper)):
             name = '%s.%s' % (self.name, name)
-            args = [name, g, getattr(self, 'run_info', None), self.source]
-            col = self.__class__(*args)
+            col = self.__class__(name=name, data=g, source=self.source,
+                                 run_info=getattr(self, 'run_info', None))
             subsets.append(col)
         return subsets
 
@@ -419,7 +419,7 @@ def _build_entity_index(self, run_info, sampling_rate):
         self.timestamps = pd.concat(_timestamps, axis=0, sort=True)
         return pd.concat(index, axis=0, sort=True).reset_index(drop=True)
 
-    def resample(self, sampling_rate, inplace=False, kind='linear'):
+    def resample(self, sampling_rate, integration_window=None, inplace=False, kind='linear'):
         '''Resample the Variable to the specified sampling rate.
 
         Parameters
@@ -436,7 +436,7 @@ def resample(self, sampling_rate, inplace=False, kind='linear'):
         '''
         if not inplace:
             var = self.clone()
-            var.resample(sampling_rate, True, kind)
+            var.resample(sampling_rate, integration_window, True, kind)
             return var
 
         if sampling_rate == self.sampling_rate:
@@ -447,18 +447,36 @@ def resample(self, sampling_rate, inplace=False, kind='linear'):
 
         self.index = self._build_entity_index(self.run_info, sampling_rate)
 
-        x = np.arange(n)
         num = len(self.index)
 
-        from scipy.interpolate import interp1d
-        f = interp1d(x, self.values.values.ravel(), kind=kind)
-        x_new = np.linspace(0, n - 1, num=num)
-        self.values = pd.DataFrame(f(x_new))
+        if integration_window is not None:
+            from scipy.sparse import lil_matrix
+            old_times = np.arange(n) / old_sr
+            new_times = np.arange(num) / sampling_rate
+            integrator = lil_matrix((num, n), dtype=np.uint8)
+            count = None
+            for i, new_time in enumerate(new_times):
+                cols = (old_times >= new_time) & (old_times < new_time + integration_window)
+                # This could be determined analytically, but dodging off-by-one errors
+                if count is None:
+                    count = np.sum(cols)
+                integrator[i, cols] = 1
+
+            old_vals = self.values.values
+            self.values = pd.DataFrame(integrator.tocsr().dot(old_vals) / count)
+        else:
+            from scipy.interpolate import interp1d
+            x = np.arange(n)
+            f = interp1d(x, self.values.values.ravel(), kind=kind)
+            x_new = np.linspace(0, n - 1, num=num)
+            self.values = pd.DataFrame(f(x_new))
+
         assert len(self.values) == len(self.index)
 
         self.sampling_rate = sampling_rate
 
-    def to_df(self, condition=True, entities=True, timing=True, sampling_rate=None):
+    def to_df(self, condition=True, entities=True, timing=True, sampling_rate=None,
+              integration_window=None):
         '''Convert to a DataFrame, with columns for name and entities.
 
         Parameters
@@ -474,7 +492,8 @@ def to_df(self, condition=True, entities=True, timing=True, sampling_rate=None):
             sampled uniformly). If False, omits them.
         '''
         if sampling_rate not in (None, self.sampling_rate):
-            return self.resample(sampling_rate).to_df(condition, entities)
+            return self.resample(sampling_rate,
+                                 integration_window=integration_window).to_df(condition, entities)
 
         df = super(DenseRunVariable, self).to_df(condition, entities)