isi distance and basic example

4 files changed, 158 insertions, 0 deletions

diff --git a/pyspike/__init__.py b/pyspike/__init__.py
new file mode 100644
index 0000000..6651eb5
--- /dev/null
+++ b/pyspike/__init__.py
@@ -0,0 +1,5 @@
+__all__ = ["function", "distances", "spikes"]
+
+from function import PieceWiseConstFunc
+from distances import isi_distance
+from spikes import spike_train_from_string
diff --git a/pyspike/distances.py b/pyspike/distances.py
new file mode 100644
index 0000000..d9790dc
--- /dev/null
+++ b/pyspike/distances.py
@@ -0,0 +1,75 @@
+""" distances.py
+
+Module containing several function to compute spike distances
+
+Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
+"""
+
+import numpy as np
+
+from pyspike import PieceWiseConstFunc
+
+def spike_train_from_string(s, sep=' '):
+    """ Converts a string of times into a SpikeTrain object.
+    Params:
+    - s: the string with (ordered) spike times
+    - sep: The separator between the time numbers.
+    Returns:
+    - array of spike times
+    """
+    return np.fromstring(s, sep=sep)
+
+def isi_distance(spikes1, spikes2, T_end, T_start=0.0):
+    """ Computes the instantaneous isi-distance S_isi (t) of the two given spike
+    trains. 
+    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).
+    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
+
+    # compute the interspike interval
+    nu1 = s1[1:]-s1[:-1]
+    nu2 = s2[1:]-s2[:-1]
+    
+    # compute the isi-distance
+    spike_events = np.empty(len(nu1)+len(nu2))
+    spike_events[0] = T_start
+    spike_events[-1] = T_end
+    # 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
+    # isi_values[0] = nu1[0]/nu2[0] - 1.0 if nu1[0] <= nu2[0] \
+    #                 else 1.0 - nu2[0]/nu1[0]
+    isi_values[0] = (nu1[0]-nu2[0])/max(nu1[0],nu2[0])
+    index1 = 0
+    index2 = 0
+    index = 1
+    while True:
+        # check which spike is next - from s1 or s2
+        if s1[index1+1] <= s2[index2+1]:
+            index1 += 1
+            if index1 >= len(nu1):
+                break
+            spike_events[index] = s1[index1]
+        else:
+            index2 += 1
+            if index2 >= len(nu2):
+                break
+            spike_events[index] = s2[index2]
+        # compute the corresponding isi-distance
+        isi_values[index] = (nu1[index1]-nu2[index2]) / \
+                            max(nu1[index1], nu2[index2])
+        index += 1
+    return PieceWiseConstFunc(spike_events, isi_values)
diff --git a/pyspike/function.py b/pyspike/function.py
new file mode 100644
index 0000000..c1de9cb
--- /dev/null
+++ b/pyspike/function.py
@@ -0,0 +1,60 @@
+""" function.py
+
+Module containing classes representing piece-wise constant and piece-wise linear
+functions.
+
+Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
+
+"""
+from __future__ import print_function
+
+import numpy as np
+
+class PieceWiseConstFunc:
+    """ A class representing a piece-wise constant function. """
+    
+    def __init__(self, x, y):
+        """ Constructs the piece-wise const function.
+        Params:
+        - x: array of length N+1 defining the edges of the intervals of the pwc
+        function.
+        - y: array of length N defining the function values at the intervals.
+        """
+        self.x = x
+        self.y = y
+
+    def get_plottable_data(self):
+        """ Returns two arrays containing x- and y-coordinates for immeditate
+        plotting of the piece-wise function.
+        """
+
+        x_plot = np.empty(2*len(self.x)-2)
+        x_plot[0] = self.x[0]
+        x_plot[1::2] = self.x[1:]
+        x_plot[2::2] = self.x[1:-1]
+        y_plot = np.empty(2*len(self.y))
+        y_plot[::2] = self.y
+        y_plot[1::2] = self.y
+
+        return x_plot, y_plot
+
+    def avrg(self):
+        """ Computes the average of the piece-wise const function:
+        a = 1/T int f(x) dx where T is the length of the interval.
+        Returns:
+        - the average a. 
+        """
+        print(self.x)
+        print(self.y)
+        return np.sum((self.x[1:]-self.x[:-1]) * self.y) / \
+            (self.x[-1]-self.x[0])
+
+    def abs_avrg(self):
+        """ Computes the average of the abs value of the piece-wise const 
+        function:
+        a = 1/T int |f(x)| dx where T is the length of the interval.
+        Returns:
+        - the average a. 
+        """
+        return np.sum((self.x[1:]-self.x[:-1]) * np.abs(self.y)) / \
+            (self.x[-1]-self.x[0])
diff --git a/pyspike/spikes.py b/pyspike/spikes.py
new file mode 100644
index 0000000..42b6501
--- /dev/null
+++ b/pyspike/spikes.py
@@ -0,0 +1,18 @@
+""" spikes.py
+
+Module containing several function to load and transform spike trains
+
+Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
+"""
+
+import numpy as np
+
+def spike_train_from_string(s, sep=' '):
+    """ Converts a string of times into a SpikeTrain object.
+    Params:
+    - s: the string with (ordered) spike times
+    - sep: The separator between the time numbers.
+    Returns:
+    - array of spike times
+    """
+    return np.fromstring(s, sep=sep)