+ add_auxiliary_spikes function incl test

3 files changed, 98 insertions, 55 deletions

diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index a5f146a..1784037 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -1,5 +1,5 @@
 __all__ = ["function", "distances", "spikes"]
 
 from function import PieceWiseConstFunc, PieceWiseLinFunc
-from distances import isi_distance, spike_distance
+from distances import add_auxiliary_spikes, isi_distance, spike_distance
 from spikes import spike_train_from_string, merge_spike_trains
diff --git a/pyspike/distances.py b/pyspike/distances.py
index a9a2cc8..10b1d3c 100644
--- a/pyspike/distances.py
+++ b/pyspike/distances.py
@@ -9,25 +9,45 @@ import numpy as np
 
 from pyspike import PieceWiseConstFunc, PieceWiseLinFunc
 
-def isi_distance(spikes1, spikes2, T_end, T_start=0.0):
+def add_auxiliary_spikes( spike_train, T_end , T_start=0.0):
+    """ Adds spikes at the beginning (T_start) and end (T_end) of the 
+    observation interval.
+    Args:
+    - spike_train: ordered array of spike times
+    - T_end: end time of the observation interval
+    - T_start: start time of the observation interval (default 0.0)
+    Returns:
+    - spike train with additional spikes at T_start and T_end.
+    """
+    assert spike_train[0] >= T_start, \
+           "Spike train has events before the given start time"
+    assert spike_train[-1] <= T_end, \
+           "Spike train has events after the given end time"
+    if spike_train[0] != T_start:
+        spike_train = np.insert(spike_train, 0, T_start)
+    if spike_train[-1] != T_end:
+        spike_train = np.append(spike_train, T_end)
+    return spike_train
+
+def isi_distance(spikes1, spikes2):
     """ Computes the instantaneous isi-distance S_isi (t) of the two given 
-    spike trains. 
+    spike trains. The spike trains are expected to have auxiliary spikes at the
+    beginning and end of the interval. Use the function add_auxiliary_spikes to
+    add those spikes to the spike train.
     Args:
-    - spikes1, spikes2: ordered arrays of spike times.
-    - T_end: end time of the observation interval.
-    - T_start: begin of the observation interval (default=0.0).
+    - spikes1, spikes2: ordered arrays of spike times with auxiliary spikes.
     Returns:
     - PieceWiseConstFunc describing the isi-distance.
     """
-    # add spikes at the beginning and end of the interval
-    s1 = np.empty(len(spikes1)+2)
-    s1[0] = T_start
-    s1[-1] = T_end
-    s1[1:-1] = spikes1
-    s2 = np.empty(len(spikes2)+2)
-    s2[0] = T_start
-    s2[-1] = T_end
-    s2[1:-1] = spikes2
+    # check for auxiliary spikes - first and last spikes should be identical
+    assert spikes1[0]==spikes2[0], \
+        "Given spike trains seems not to have auxiliary spikes!"
+    assert spikes1[-1]==spikes2[-1], \
+        "Given spike trains seems not to have auxiliary spikes!"
+
+    # shorter names
+    s1 = spikes1
+    s2 = spikes2
 
     # compute the interspike interval
     nu1 = s1[1:]-s1[:-1]
@@ -35,7 +55,7 @@ def isi_distance(spikes1, spikes2, T_end, T_start=0.0):
     
     # compute the isi-distance
     spike_events = np.empty(len(nu1)+len(nu2))
-    spike_events[0] = T_start
+    spike_events[0] = s1[0]
     # the values have one entry less - the number of intervals between events
     isi_values = np.empty(len(spike_events)-1)
     # add the distance of the first events
@@ -69,7 +89,7 @@ def isi_distance(spikes1, spikes2, T_end, T_start=0.0):
                             max(nu1[index1], nu2[index2])
         index += 1
     # the last event is the interval end
-    spike_events[index] = T_end
+    spike_events[index] = s1[-1]
     # use only the data added above 
     # could be less than original length due to equal spike times
     return PieceWiseConstFunc(spike_events[:index+1], isi_values[:index])
@@ -77,7 +97,7 @@ def isi_distance(spikes1, spikes2, T_end, T_start=0.0):
 
 def get_min_dist(spike_time, spike_train, start_index=0):
     """ Returns the minimal distance |spike_time - spike_train[i]| 
-    with i>=start_index
+    with i>=start_index.
     """
     d = abs(spike_time - spike_train[start_index])
     start_index += 1
@@ -91,31 +111,28 @@ def get_min_dist(spike_time, spike_train, start_index=0):
     return d
 
 
-def spike_distance(spikes1, spikes2, T_end, T_start=0.0):
+def spike_distance(spikes1, spikes2):
     """ Computes the instantaneous spike-distance S_spike (t) of the two given
-    spike trains.
+    spike trains. The spike trains are expected to have auxiliary spikes at the
+    beginning and end of the interval. Use the function add_auxiliary_spikes to
+    add those spikes to the spike train.
     Args:
-    - spikes1, spikes2: ordered arrays of spike times.
-    - T_end: end time of the observation interval.
-    - T_start: begin of the observation interval (default=0.0).
+    - spikes1, spikes2: ordered arrays of spike times with auxiliary spikes.
     Returns:
     - PieceWiseLinFunc describing the spike-distance.
     """
-    # add spikes at the beginning and end of the interval
-    t1 = np.empty(len(spikes1)+2)
-    t1[0] = T_start
-    t1[-1] = T_end
-    t1[1:-1] = spikes1
-    t2 = np.empty(len(spikes2)+2)
-    t2[0] = T_start
-    t2[-1] = T_end
-    t2[1:-1] = spikes2
+    # check for auxiliary spikes - first and last spikes should be identical
+    assert spikes1[0]==spikes2[0], \
+        "Given spike trains seems not to have auxiliary spikes!"
+    assert spikes1[-1]==spikes2[-1], \
+        "Given spike trains seems not to have auxiliary spikes!"
+    # shorter variables
+    t1 = spikes1
+    t2 = spikes2
 
     spike_events = np.empty(len(t1)+len(t2)-2)
-    spike_events[0] = T_start
-    spike_events[-1] = T_end
+    spike_events[0] = t1[0]
     y_starts = np.empty(len(spike_events)-1)
-    y_starts[0] = 0.0
     y_ends = np.empty(len(spike_events)-1)
 
     index1 = 0
@@ -125,10 +142,13 @@ def spike_distance(spikes1, spikes2, T_end, T_start=0.0):
     dt_f1 = get_min_dist(t1[1], t2, 0)
     dt_p2 = 0.0
     dt_f2 = get_min_dist(t2[1], t1, 0)
-    isi1 = t1[1]-t1[0]
-    isi2 = t2[1]-t2[0]
+    isi1 = max(t1[1]-t1[0], t1[2]-t1[1])
+    isi2 = max(t2[1]-t2[0], t2[2]-t2[1])
+    s1 = dt_f1*(t1[1]-t1[0])/isi1
+    s2 = dt_f2*(t2[1]-t2[0])/isi2
+    y_starts[0] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
     while True:
-        print(index, index1, index2)
+        # print(index, index1, index2)
         if t1[index1+1] < t2[index2+1]:
             index1 += 1
             # break condition relies on existence of spikes at T_end
@@ -166,8 +186,8 @@ def spike_distance(spikes1, spikes2, T_end, T_start=0.0):
             index2 += 1
             if (index1+1 >= len(t1)) or (index2+1 >= len(t2)):
                 break
-            assert( dt_f2 == 0.0 )
-            assert( dt_f1 == 0.0 )
+            assert dt_f2 == 0.0
+            assert dt_f1 == 0.0
             spike_events[index] = t1[index1]
             y_ends[index-1] = 0.0
             y_starts[index] = 0.0
@@ -179,9 +199,13 @@ def spike_distance(spikes1, spikes2, T_end, T_start=0.0):
             isi2 = t2[index2+1]-t2[index2]
         index += 1
     # the last event is the interval end
-    spike_events[index] = T_end
-    # the ending value of the last interval is 0
-    y_ends[index-1] = 0.0
+    spike_events[index] = t1[-1]
+    # the ending value of the last interval
+    isi1 = max(t1[-1]-t1[-2], t1[-2]-t1[-3])
+    isi2 = max(t2[-1]-t2[-2], t2[-2]-t2[-3])
+    s1 = dt_p1*(t1[-1]-t1[-2])/isi1
+    s2 = dt_p2*(t2[-1]-t2[-2])/isi2
+    y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
     # use only the data added above 
     # could be less than original length due to equal spike times
     return PieceWiseLinFunc(spike_events[:index+1], 
diff --git a/test/test_distance.py b/test/test_distance.py
index 17ca14a..35bdf85 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -11,6 +11,13 @@ from numpy.testing import assert_equal, assert_array_almost_equal
 
 import pyspike as spk
 
+def test_auxiliary_spikes():
+    t = np.array([0.2, 0.4, 0.6, 0.7])
+    t_aux = spk.add_auxiliary_spikes(t, T_end=1.0, T_start=0.1)
+    assert_equal(t_aux, [0.1, 0.2, 0.4, 0.6, 0.7, 1.0])
+    t_aux = spk.add_auxiliary_spikes(t_aux, 1.0)
+    assert_equal(t_aux, [0.0, 0.1, 0.2, 0.4, 0.6, 0.7, 1.0])
+
 def test_isi():
     # generate two spike trains:
     t1 = np.array([0.2, 0.4, 0.6, 0.7])
@@ -21,7 +28,11 @@ def test_isi():
     expected_isi = [-0.1/0.3, -0.1/0.3, 0.05/0.2, 0.05/0.2, -0.15/0.35, 
                     -0.25/0.35, -0.05/0.35, 0.2/0.3, 0.25/0.3, 0.25/0.3]
     
-    f = spk.isi_distance(t1, t2, 1.0)
+    t1 = spk.add_auxiliary_spikes(t1, 1.0)
+    t2 = spk.add_auxiliary_spikes(t2, 1.0)
+    f = spk.isi_distance(t1, t2)
+
+    print("ISI: ", f.y)
 
     assert_equal(f.x, expected_times)
     assert_array_almost_equal(f.y, expected_isi, decimal=14)
@@ -33,7 +44,9 @@ def test_isi():
     expected_times = [0.0,0.1,0.2,0.4,0.5,0.6,1.0]
     expected_isi = [0.1/0.2, -0.1/0.3, -0.1/0.3, 0.1/0.2, 0.1/0.2, -0.0/0.5]
 
-    f = spk.isi_distance(t1, t2, 1.0)
+    t1 = spk.add_auxiliary_spikes(t1, 1.0)
+    t2 = spk.add_auxiliary_spikes(t2, 1.0)
+    f = spk.isi_distance(t1, t2)
 
     assert_equal(f.x, expected_times)
     assert_array_almost_equal(f.y, expected_isi, decimal=14)
@@ -46,16 +59,19 @@ def test_spike():
 
     # pen&paper calculation of the spike distance
     expected_times = [0.0,0.2,0.3,0.4,0.45,0.6,0.7,0.8,0.9,0.95,1.0]
-    s1 = np.array([0.0, 0.1, (0.1*0.1+0.05*0.1)/0.2, 0.05, (0.05*0.15 * 2)/0.2,
-          0.15, 0.1, 0.1*0.2/0.3, 0.1**2/0.3, 0.1*0.05/0.3, 0.0])
-    s2 = np.array([0.0, 0.1*0.2/0.3, 0.1, (0.1*0.05 * 2)/.15, 0.05, 
-          (0.05*0.2+0.1*0.15)/0.35, (0.05*0.1+0.1*0.25)/0.35, 0.1,0.1,0.05,0.0])
+    s1 = np.array([0.1, 0.1, (0.1*0.1+0.05*0.1)/0.2, 0.05, (0.05*0.15 * 2)/0.2,
+                   0.15, 0.1, 0.1*0.2/0.3, 0.1**2/0.3, 0.1*0.05/0.3, 0.1])
+    s2 = np.array([0.1, 0.1*0.2/0.3, 0.1, (0.1*0.05 * 2)/.15, 0.05, 
+                   (0.05*0.2+0.1*0.15)/0.35, (0.05*0.1+0.1*0.25)/0.35, 
+                   0.1, 0.1, 0.05, 0.05])
     isi1 = np.array([0.2, 0.2, 0.2, 0.2, 0.2, 0.1, 0.3, 0.3, 0.3, 0.3])
     isi2 = np.array([0.3, 0.3, 0.15, 0.15, 0.35, 0.35, 0.35, 0.1, 0.05, 0.05])
     expected_y1 = (s1[:-1]*isi2+s2[:-1]*isi1) / (0.5*(isi1+isi2)**2)
     expected_y2 = (s1[1:]*isi2+s2[1:]*isi1) / (0.5*(isi1+isi2)**2)
 
-    f = spk.spike_distance(t1, t2, 1.0)
+    t1 = spk.add_auxiliary_spikes(t1, 1.0)
+    t2 = spk.add_auxiliary_spikes(t2, 1.0)
+    f = spk.spike_distance(t1, t2)
 
     assert_equal(f.x, expected_times)
     assert_array_almost_equal(f.y1, expected_y1, decimal=14)
@@ -66,20 +82,23 @@ def test_spike():
     t2 = np.array([0.1,0.4,0.5,0.6])
 
     expected_times = [0.0,0.1,0.2,0.4,0.5,0.6,1.0]
-    s1 = np.array([0.0, 0.1*0.1/0.2, 0.1, 0.0, 0.0, 0.0, 0.0])
-    s2 = np.array([0.0, 0.1, 0.1*0.2/0.3, 0.0, 0.1, 0.0, 0.0])
+    s1 = np.array([0.1, 0.1*0.1/0.2, 0.1, 0.0, 0.0, 0.0, 0.0])
+    s2 = np.array([0.1*0.1/0.3, 0.1, 0.1*0.2/0.3, 0.0, 0.1, 0.0, 0.0])
     isi1 = np.array([0.2, 0.2, 0.2, 0.2, 0.2, 0.4])
-    isi2 = np.array([0.1, 0.3, 0.3, 0.1, 0.1, 0.4])
+    isi2 = np.array([0.3, 0.3, 0.3, 0.1, 0.1, 0.4])
     expected_y1 = (s1[:-1]*isi2+s2[:-1]*isi1) / (0.5*(isi1+isi2)**2)
     expected_y2 = (s1[1:]*isi2+s2[1:]*isi1) / (0.5*(isi1+isi2)**2)
     
-    f = spk.spike_distance(t1, t2, 1.0)
+    t1 = spk.add_auxiliary_spikes(t1, 1.0)
+    t2 = spk.add_auxiliary_spikes(t2, 1.0)
+    f = spk.spike_distance(t1, t2)
 
     assert_equal(f.x, expected_times)
     assert_array_almost_equal(f.y1, expected_y1, decimal=14)
     assert_array_almost_equal(f.y2, expected_y2, decimal=14)
 
 
-if __name__ == "main":
+if __name__ == "__main__":
+    test_auxiliary_spikes()
     test_isi()
     test_spike()