Merge pull request #571 from effigies/natural_density

ENH: Calculate default sampling rate for SparseRunVariable.to_dense

bids/variables/tests/test_variables.py

@@ -82,6 +82,29 @@ def test_sparse_run_variable_to_dense(layout1):
         assert dense.source == 'events'
 
 
+def test_sparse_run_variable_to_dense_default_sr(layout1):
+    index = load_variables(layout1, types='events', scan_length=480)
+    runs = index.get_nodes('run', {'subject': ['01', '02']})
+
+    for i, run in enumerate(runs):
+        var = run.variables['RT']
+        dense = var.to_dense()
+
+        # Check that a sensible sampling rate was found
+        assert np.allclose(dense.sampling_rate, 1)
+
+        # Check that all unique values are identical
+        sparse_vals = set(np.unique(var.values.values)) | {0}
+        dense_vals = set(np.unique(dense.values.values))
+        assert sparse_vals == dense_vals
+
+        assert len(dense.values) > len(var.values)
+        assert isinstance(dense, DenseRunVariable)
+        assert dense.values.shape == (480, 1)
+        assert len(dense.run_info) == len(var.run_info)
+        assert dense.source == 'events'
+
+
 def test_merge_densified_variables(layout1):
     SR = 10
     dataset = load_variables(layout1, types='events', scan_length=480)

bids/variables/variables.py

@@ -8,6 +8,7 @@
 from abc import abstractmethod, ABCMeta
 from bids.utils import listify
 from itertools import chain
+from functools import reduce
 from bids.utils import matches_entities
 
 class BIDSVariable(metaclass=ABCMeta):
@@ -362,19 +363,28 @@ def get_duration(self):
         """Return the total duration of the Variable's run(s). """
         return sum([r.duration for r in self.run_info])
 
-    def to_dense(self, sampling_rate):
+    def to_dense(self, sampling_rate=None):
         """Convert the current sparse column to a dense representation.
 
+        If sampling_rate is not provided, the largest interval able to
+        faithfully represent all onsets and durations will be determined.
+        The sampling rate is the reciprocal of that interval.
+
         Parameters
         ----------
-        sampling_rate : int or str
-            Sampling rate (in Hz) to use when
-            constructing the DenseRunVariable.
+        sampling_rate : float or None
+            Sampling rate (in Hz) to use when constructing the DenseRunVariable.
 
         Returns
         -------
         DenseRunVariable
         """
+        if sampling_rate is None:
+            # Cast onsets and durations to milliseconds
+            onsets = np.round(self.onset * 1000).astype(int)
+            durations = np.round(self.duration * 1000).astype(int)
+            sampling_rate = 1000. / reduce(math.gcd, [*onsets, *durations])
+
         duration = int(math.ceil(sampling_rate * self.get_duration()))
         ts = np.zeros(duration, dtype=self.values.dtype)