[WIP] Enh/58 interpolation methods

neuroglia/event.py

@@ -1,13 +1,19 @@
-import numpy as np
 import pandas as pd
 import xarray as xr
+from six import string_types
+from collections import namedtuple
 
-from sklearn.base import BaseEstimator,TransformerMixin
+from sklearn.base import BaseEstimator, TransformerMixin
+from scipy.interpolate import InterpolatedUnivariateSpline
 
-from .utils import create_interpolator, events_to_xr_dim
+from .utils import events_to_xr_dim
 from .spike import Binner, DEFAULT_TAU
+from .interpolator import kriging_interpolator_factory, \
+    sinc_interpolator_factory
+
+
+class PeriEventTraceSampler(BaseEstimator, TransformerMixin):
 
-class PeriEventTraceSampler(BaseEstimator,TransformerMixin):
     """Take event-aligned samples of traces from a population of neurons.
 
     Traces are sampled relative to the event time. There is no enforced
@@ -17,24 +23,55 @@ class PeriEventTraceSampler(BaseEstimator,TransformerMixin):
 
     Parameters
     ----------
-    traces : pandas DataFrame with 'time' as the index and neuron IDs in columns
-        The traces that will be sampled from when the transform method is called
+    traces : pandas DataFrame with 'time' as the index and neuron IDs in
+    columns
+        The traces that will be sampled from when the transform method is
+        called
     sample_times : array
         Time relative to events that will be used to sample or bin spikes.
+    interpolator : string or callable, default: `cubic_spline`
+        interpolator factory or string correponding to a builtin interpolator;
+        builtin interpolators are: "cubic_spline", "kriging", "sinc"
 
     Notes
     -----
-
-    This estimator is stateless (besides constructor parameters), the
+    - This estimator is stateless (besides constructor parameters), the
     fit method does nothing but is useful when used in a pipeline.
     """
-    def __init__(self, traces, sample_times):
+
+    Interpolators = namedtuple(
+        "BuiltinInterpolators",
+        ["cubic_spline", "kriging", "sinc", ]
+    )(
+        InterpolatedUnivariateSpline,
+        kriging_interpolator_factory,
+        sinc_interpolator_factory,
+    )  # TODO figure out a better name...
+
+    def __init__(
+            self,
+            traces,
+            sample_times,
+            interpolator=Interpolators.cubic_spline,
+    ):
         self.sample_times = sample_times
         self.traces = traces
 
+        if isinstance(interpolator, string_types):
+            try:
+                self.interpolator = getattr(
+                    PeriEventSpikeSampler.Interpolators
+                )
+            except AttributeError:
+                raise ValueError(
+                    "unknown builtin interpolator: " + interpolator
+                )
+        else:
+            self.interpolator = interpolator
+
     def _make_splined_traces(self):
         self.splined_traces_ = self.traces.apply(
-            lambda y: create_interpolator(self.traces.index,y),
+            lambda y: self.interpolator(self.traces.index, y),
             axis=0,
         )
 
@@ -72,16 +109,19 @@ def transform(self, X):
         def extractor(ev):
             t = self.sample_times + ev['time']
             interpolated = self.splined_traces_.apply(
-                lambda s: pd.Series(s(t),index=self.sample_times)
-                )
-            return xr.DataArray(interpolated.T,dims=['sample_times','neuron'])
+                lambda s: pd.Series(s(t), index=self.sample_times)
+            )
+            return xr.DataArray(
+                interpolated.T,
+                dims=['sample_times', 'neuron', ]
+            )
 
         # do the extraction
-        tensor = [extractor(ev) for _,ev in X.iterrows()]
+        tensor = [extractor(ev) for _, ev in X.iterrows()]
         concat_dim = events_to_xr_dim(X)
 
         # concatenate the DataArrays into a single DataArray
-        return xr.concat(tensor,dim=concat_dim)
+        return xr.concat(tensor, dim=concat_dim)
 
 
 class PeriEventSpikeSampler(BaseEstimator,TransformerMixin):

neuroglia/interpolator.py

@@ -0,0 +1,86 @@
+import numpy as np
+from sklearn.gaussian_process import GaussianProcessRegressor
+
+
+def _sinc_interpolator(x, s, u):
+    """sinc interpolator
+
+    from:
+    https://gist.github.com/endolith/1297227
+
+    Parameters
+    ----------
+    x : `np.ndarray`
+        sampled values
+    s : `np.ndarray`
+        sampled points
+    u : `np.ndarray`
+        unsampled points
+
+    Returns
+    -------
+    np.ndarray
+        interpolated values
+
+    Notes
+    -----
+    - `x` is typically the dependent variable and `s` is typically the
+    independent variable
+    """
+    if len(x) != len(s):
+        raise ValueError('x and s must be the same length')
+
+    # Find the period
+    T = s[1] - s[0]
+
+    sincM = np.tile(u, (len(s), 1)) - np.tile(s[:, np.newaxis], (1, len(u)))
+
+    return np.dot(x, np.sinc(sincM / T))
+
+
+def kriging_interpolator_factory(x, y, **kwargs):
+    """1-dimensional kriging
+
+    from:
+    https://stackoverflow.com/questions/24978052/interpolation-over-regular-grid-in-python
+
+    Parameters
+    ----------
+    x : `np.ndarray`
+        independent variable used to generate interpolator
+    y : `np.ndarray`
+        dependent variable used to generator interpolator
+    **kwargs
+        keyword arguments supplied to `GaussianProcessRegressor`
+
+    Returns
+    -------
+    function
+        interpolator function
+
+    Notes
+    -----
+    - assumes `x` and `y` are of shape: (n, ), with an equal number of elements
+    """
+    process = GaussianProcessRegressor(**kwargs)
+    process.fit(x.reshape(1, -1), y.reshape(1, -1))
+
+    return lambda u: process.predict(u.reshape(1, -1)).reshape(-1)
+
+
+def sinc_interpolator_factory(x, y):
+    """sinc interpolator factory
+
+    Parameters
+    ----------
+    x : `np.ndarray`
+        independent variable used to generate interpolator
+    y : `np.ndarray`
+        dependent variable used to generator interpolator
+
+    Returns
+    -------
+    function
+        interpolator function
+    """
+    return lambda u: _sinc_interpolator(y, x, u)

neuroglia/utils.py

@@ -1,30 +1,9 @@
-import numpy as np
 import xarray as xr
-from scipy import interpolate
 
-def events_to_xr_dim(events,dim='event'):
+
+def events_to_xr_dim(events, dim='event'):
     # builds the event dataframe into coords for xarray
     coords = events.to_dict(orient='list')
-    coords = {k:(dim,v) for k,v in coords.items()}
+    coords = {k: (dim, v, ) for k, v in coords.items()}
     # define a DataArray that will describe the event dimension
-    return xr.DataArray(
-        events.index,
-        name=dim,
-        dims=[dim],
-        coords=coords,
-    )
-
-def create_interpolator(t,y):
-    """ creates a cubic spline interpolator
-
-    Parameters
-    ---------
-    y : array-like of floats
-
-    Returns
-    ---------
-    interpolator function that accepts a list of times
-
-    """
-    interpolator = interpolate.InterpolatedUnivariateSpline(t, y)
-    return interpolator
+    return xr.DataArray(events.index, name=dim, dims=[dim], coords=coords)

requirements.txt

@@ -3,3 +3,4 @@ scipy
 pandas
 xarray
 scikit-learn
+six

tests/test_interpolator.py

@@ -0,0 +1,64 @@
+import pytest
+import numpy as np
+import numpy.testing as npt
+
+from neuroglia import interpolator
+
+
+# values from justin kiggins
+TIME_PRECISE = np.arange(0, 0.2, 0.001)
+TIME_TARGET = TIME_PRECISE[::10]  # real time
+TIME_OBSERVED = TIME_TARGET + 0.002 * np.random.randn(*TIME_TARGET.shape)  # unevenly sampled time
+make_neuron = lambda t: np.sin(100 * t)
+
+
+@pytest.mark.parametrize("x, y, u, expected", [
+    (
+        np.linspace(-np.pi, np.pi, 201),
+        np.sin(np.linspace(-np.pi, np.pi, 201)),
+        np.linspace(-np.pi, np.pi, 201),
+        np.sin(np.linspace(-np.pi, np.pi, 201)),
+    ),
+    (
+        np.linspace(0, 5, 100),
+        np.linspace(0, 5, 100),
+        np.linspace(0, 5, 100),
+        np.linspace(0, 5, 100),
+    ),
+    (
+        TIME_OBSERVED,
+        make_neuron(TIME_OBSERVED),
+        TIME_OBSERVED,
+        make_neuron(TIME_OBSERVED),
+    ),
+])
+def test_kriging_interpolator_factory(x, y, u, expected):
+    """Bad test...
+    """
+    npt.assert_allclose(
+        interpolator.kriging_interpolator_factory(x, y)(u),
+        expected,
+        atol=0.05
+    )
+
+
+@pytest.mark.parametrize("x, y, u, expected", [
+    (
+        np.linspace(-np.pi, np.pi, 201),
+        np.sin(np.linspace(-np.pi, np.pi, 201)),
+        np.linspace(-np.pi, np.pi, 201),
+        np.sin(np.linspace(-np.pi, np.pi, 201)),
+    ),
+    (
+        np.linspace(0, 5, 100),
+        np.linspace(0, 5, 100),
+        np.linspace(0, 5, 100),
+        np.linspace(0, 5, 100),
+    ),
+])
+def test_sinc_interpolator_factory(x, y, u, expected):
+    npt.assert_allclose(
+        interpolator.sinc_interpolator_factory(x, y)(u),
+        expected,
+        atol=0.05
+    )