Repository: bids-standard/pybids

Issue:
Lag transformation
Thinking of implementing Lag, since it's such a common transformation in practice.

----

Current docs:

> ### Lag(Input, Output=None):
> [TODO]
> 
> Arguments:
> * Input(string, mandatory): The name of the variable to operate on.
> * Output (string, optional): the optional list of column names to write out to. By default, computation is done in-place (i.e., input columnise overwritten).

----

It's a bit under-specified, so while fleshing it out (and implementing) we have opportunities to think about what would make it most useful. Questions:

1) Should we permit a lag != 1? If so, is definitially `Lag(Lag(x, 1), 1) = Lag(x, 2)`?
2) Should we permit negative lags?
3) Should we allow taking the difference in the same step? So we get `Lag(x, diff=True) = Lag(x, diff=False) - x`? (or `x - Lag(x, diff=False)`)?
4) Should we default to overwriting columns or generating new columns from the original name? BEP 012 [would suggest](https://github.com/bids-standard/bids-specification/pull/519/files#diff-093ab536060c8ee949aa452200ec2676377a2fcd1a8674f91c8e6a1cfd42281fR127) calling `Lag(x) = x_shift_back`, though `x_lag` might be clearer.

Any other thoughts?

cc @adelavega 

Discussion:

Patch:
diff --git a/bids/modeling/transformations/__init__.py b/bids/modeling/transformations/__init__.py
index c2f03356..4d642384 100644
--- a/bids/modeling/transformations/__init__.py
+++ b/bids/modeling/transformations/__init__.py
@@ -1,5 +1,5 @@
 from .compute import (Sum, Product, Scale, Orthogonalize, Threshold, And_, Or_,
-                      Not, Demean, Convolve)
+                      Not, Demean, Convolve, Lag)
 from .munge import (Split, Rename, Assign, Copy, Factor, Filter, Select,
                     Delete, DropNA, Replace, ToDense, Group, Resample)
 from .base import TransformerManager
@@ -19,6 +19,7 @@ __all__ = [
     'Factor',
     'Filter',
     'Group',
+    'Lag',
     'Not',
     'Or',
     'Or_',
diff --git a/bids/modeling/transformations/compute.py b/bids/modeling/transformations/compute.py
index 8a6cba8f..9478f926 100644
--- a/bids/modeling/transformations/compute.py
+++ b/bids/modeling/transformations/compute.py
@@ -4,6 +4,7 @@ Transformations that primarily involve numerical computation on variables.
 import math
 import numpy as np
 import pandas as pd
+from scipy import ndimage as ndi
 from bids.utils import listify
 from .base import Transformation
 from bids.modeling import hrf
@@ -290,3 +291,68 @@ class Or_(Transformation):
     def _transform(self, dfs):
         df = pd.concat(dfs, axis=1, sort=True)
         return df.any(axis=1).astype(int)
+
+
+class Lag(Transformation):
+    """Lag variable
+
+    Returns a variable that is lagged by a specified number of time points.
+    Spline interpolation of the requested ``order`` is used for non-integer
+    shifts.
+    Points outside the input are filled according to the given ``mode``.
+    Negative shifts move values toward the beginning of the sequence.
+    If ``difference`` is ``True``, the backward difference is calculated for
+    positive shifts, and the forward difference is calculated for negative
+    shifts.
+
+    Additional ``mode``s may be defined if there is need for them.
+    The ``scipy.ndimage.shift`` method provides the reference implementation
+    for all current modes. "Constant" is equivalent to the shift parameter
+    "cval".
+
+
+    Parameters
+    ----------
+    var : :obj:`numpy.ndarray`
+        variable to lag
+    shift : float, optional
+        number of places to shift the values (default: 1)
+    order : int, optional
+        order of spline interpolation, from 0-5 (default: 3)
+    mode : string
+        the `mode` parameter determines how the input array is extended
+        beyond its boundaries. Default is 'nearest'.
+        The following values are accepted:
+
+        "nearest" (a a a a | a b c d | d d d d)
+            The input is extended by replicating the boundary values
+        "reflect" (d c b a | a b c d | d c b a)
+            The input is extended by reflecting the array over the edge
+        "constant" (k k k k | a b c d | k k k k)
+            The input is extended by filling all values beyond the edge
+            with the same constant value, defined by the ``constant`` parameter
+    constant : float, optional
+        value to fill past edges of input if ``mode`` is ``"constant"``.
+        (default: 0)
+    difference : boolean, optional
+        Calculate the backward (if shift is positive) or forward (if shift is
+        negative) difference.
+        For the forward difference dx of an array x, dx[i] = x[i+1] - x[i].
+        For the backward difference dx of an array x, dx[i] = x[i] - x[i-1].
+        (default: ``False``)
+    """
+
+    _input_type = 'numpy'
+    _return_type = 'numpy'
+
+    def _transform(self, var, shift=1, order=3, mode="nearest",
+                   constant=0.0, difference=False):
+        var = var.flatten()
+        shifted = ndi.shift(var, shift=shift, order=order, mode=mode,
+                            cval=constant)
+        if not difference:
+            return shifted
+        elif shift >= 0:
+            return var - shifted
+        else:
+            return shifted - var


Test patch:
diff --git a/bids/modeling/tests/test_transformations.py b/bids/modeling/tests/test_transformations.py
index 8cdcf84b..2963de54 100644
--- a/bids/modeling/tests/test_transformations.py
+++ b/bids/modeling/tests/test_transformations.py
@@ -1,6 +1,6 @@
 # from bids.modeling.variables import load_variables
 from bids.modeling import transformations as transform
-from bids.variables import SparseRunVariable
+from bids.variables import SparseRunVariable, DenseRunVariable
 from bids.variables.entities import RunInfo
 from bids.variables.collections import BIDSRunVariableCollection
 from bids.layout import BIDSLayout
@@ -542,3 +542,56 @@ def test_resample(collection):
     # Assert that the auc is more or less the same (not exact, rounding error)
     # Values are around 0.25
     assert np.allclose(old_auc, new_auc, rtol=0.05)
+
+
+def test_Lag():
+    var = DenseRunVariable(
+        name="rot_x",
+        values=np.arange(5., 20.),
+        run_info=RunInfo({}, 15, 1, "none", 15),
+        source='regressors',
+        sampling_rate=1
+    )
+    coll = BIDSRunVariableCollection([var], sampling_rate=1)
+
+    # Forward shift
+    transform.Lag(coll, "rot_x", output="d_rot_x")
+    d_rot_x = coll["d_rot_x"].values.values
+    assert np.isclose(d_rot_x[0, 0], 5.)
+    assert np.allclose(d_rot_x[1:, 0], np.arange(5., 19.))
+
+    # Backward shift
+    transform.Lag(coll, "rot_x", output="d_rot_x", shift=-1)
+    d_rot_x = coll["d_rot_x"].values.values
+    assert np.isclose(d_rot_x[-1, 0], 19.)
+    assert np.allclose(d_rot_x[:-1, 0], np.arange(6., 20.))
+
+    # Half shift; don't know why you'd want to do it, but you can
+    transform.Lag(coll, "rot_x", output="half_shift", shift=0.5, order=1)
+    half_shift = coll["half_shift"].values.values
+    assert np.isclose(half_shift[0, 0], 5.)
+    assert np.allclose(half_shift[1:, 0], np.arange(5.5, 19.5))
+
+    # Constant mode
+    transform.Lag(coll, "rot_x", output="d_rot_x", mode="constant")
+    d_rot_x = coll["d_rot_x"].values.values
+    assert np.isclose(d_rot_x[0, 0], 0.)
+    assert np.allclose(d_rot_x[1:, 0], np.arange(5., 19.))
+
+    # Reflect mode
+    transform.Lag(coll, "rot_x", output="d_rot_x", mode="reflect")
+    d_rot_x = coll["d_rot_x"].values.values
+    assert np.isclose(d_rot_x[0, 0], 5.)
+    assert np.allclose(d_rot_x[1:, 0], np.arange(5., 19.))
+
+    # Forward shift -> Backward difference
+    transform.Lag(coll, "rot_x", output="d_rot_x", difference=True)
+    d_rot_x = coll["d_rot_x"].values.values
+    assert np.isclose(d_rot_x[0, 0], 0.)
+    assert np.allclose(d_rot_x[1:, 0], 1.)
+
+    # Backward shift -> Forward difference
+    transform.Lag(coll, "rot_x", output="d_rot_x", shift=-1, difference=True)
+    d_rot_x = coll["d_rot_x"].values.values
+    assert np.isclose(d_rot_x[-1, 0], 0.)
+    assert np.allclose(d_rot_x[:-1, 0], 1.)