+average_profile function

3 files changed, 78 insertions, 1 deletions

diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index c58a6b1..1bfa7fc 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -6,7 +6,7 @@ Distributed under the BSD License
 
 __all__ = ["function", "distances", "spikes"]
 
-from function import PieceWiseConstFunc, PieceWiseLinFunc
+from function import PieceWiseConstFunc, PieceWiseLinFunc, average_profile
 from distances import isi_distance, spike_distance, \
     isi_distance_multi, spike_distance_multi, isi_distance_matrix
 from spikes import add_auxiliary_spikes, load_spike_trains_from_txt, \
diff --git a/pyspike/function.py b/pyspike/function.py
index bd3e2d5..46fdea2 100644
--- a/pyspike/function.py
+++ b/pyspike/function.py
@@ -26,9 +26,17 @@ class PieceWiseConstFunc:
         function.
         - y: array of length N defining the function values at the intervals.
         """
+        # convert parameters to arrays, also ensures copying
         self.x = np.array(x)
         self.y = np.array(y)
 
+    def copy(self):
+        """ Returns a copy of itself
+        Returns:
+        - PieceWiseConstFunc copy
+        """
+        return PieceWiseConstFunc(self.x, self.y)
+
     def almost_equal(self, other, decimal=14):
         """ Checks if the function is equal to another function up to `decimal`
         precision.
@@ -108,10 +116,18 @@ class PieceWiseLinFunc:
         - y2: array of length N defining the function values at the right of
         the intervals.
         """
+        # convert to array, which also ensures copying
         self.x = np.array(x)
         self.y1 = np.array(y1)
         self.y2 = np.array(y2)
 
+    def copy(self):
+        """ Returns a copy of itself
+        Returns:
+        - PieceWiseLinFunc copy
+        """
+        return PieceWiseLinFunc(self.x, self.y1, self.y2)
+
     def almost_equal(self, other, decimal=14):
         """ Checks if the function is equal to another function up to `decimal`
         precision.
@@ -183,3 +199,22 @@ class PieceWiseLinFunc:
         """
         self.y1 *= fac
         self.y2 *= fac
+
+
+def average_profile(profiles):
+    """ Computes the average profile from the given ISI- or SPIKE-profiles.
+    Args:
+    - profiles: list of PieceWiseConstFunc or PieceWiseLinFunc representing
+    ISI- or SPIKE-profiles to be averaged
+    Returns:
+    - avrg_profile: PieceWiseConstFunc or PieceWiseLinFunc containing the
+    average profile.
+    """
+    assert len(profiles) > 1
+
+    avrg_profile = profiles[0].copy()
+    for i in xrange(1, len(profiles)):
+        avrg_profile.add(profiles[i])
+    avrg_profile.mul_scalar(1.0/len(profiles))  # normalize
+
+    return avrg_profile
diff --git a/test/test_function.py b/test/test_function.py
index c5caa5a..ed70180 100644
--- a/test/test_function.py
+++ b/test/test_function.py
@@ -68,6 +68,23 @@ def test_pwc_mul():
     assert_array_almost_equal(f.y, 1.5/5.0*np.array(y), decimal=16)
 
 
+def test_pwc_avrg():
+    # some random data
+    x = [0.0, 1.0, 2.0, 2.5, 4.0]
+    y = [1.0, -0.5, 1.5, 0.75]
+    f1 = spk.PieceWiseConstFunc(x, y)
+
+    x = [0.0, 0.75, 2.0, 2.5, 2.7, 4.0]
+    y = [0.5, 1.0, -0.25, 0.0, 1.5]
+    f2 = spk.PieceWiseConstFunc(x, y)
+
+    f_avrg = spk.average_profile([f1, f2])
+    x_expected = [0.0, 0.75, 1.0, 2.0, 2.5, 2.7, 4.0]
+    y_expected = [0.75, 1.0, 0.25, 0.625, 0.375, 1.125]
+    assert_array_almost_equal(f_avrg.x, x_expected, decimal=16)
+    assert_array_almost_equal(f_avrg.y, y_expected, decimal=16)
+
+
 def test_pwl():
     x = [0.0, 1.0, 2.0, 2.5, 4.0]
     y1 = [1.0, -0.5, 1.5, 0.75]
@@ -134,6 +151,31 @@ def test_pwl_mul():
     assert_array_almost_equal(f.y1, 1.5/5.0*np.array(y1), decimal=16)
     assert_array_almost_equal(f.y2, 1.5/5.0*np.array(y2), decimal=16)
 
+
+def test_pwl_avrg():
+    x = [0.0, 1.0, 2.0, 2.5, 4.0]
+    y1 = [1.0, -0.5, 1.5, 0.75]
+    y2 = [1.5, -0.4, 1.5, 0.25]
+    f1 = spk.PieceWiseLinFunc(x, y1, y2)
+
+    x = [0.0, 0.75, 2.0, 2.5, 2.7, 4.0]
+    y1 = [0.5, 1.0, -0.25, 0.0, 1.5]
+    y2 = [0.8, 0.2, -1.0, 0.0, 2.0]
+    f2 = spk.PieceWiseLinFunc(x, y1, y2)
+
+    x_expected = [0.0, 0.75, 1.0, 2.0, 2.5, 2.7, 4.0]
+    y1_expected = np.array([1.5, 1.0+1.0+0.5*0.75, -0.5+1.0-0.8*0.25/1.25,
+                            1.5-0.25, 0.75, 1.5+0.75-0.5*0.2/1.5]) / 2
+    y2_expected = np.array([0.8+1.0+0.5*0.75, 1.5+1.0-0.8*0.25/1.25, -0.4+0.2,
+                            1.5-1.0, 0.75-0.5*0.2/1.5, 2.25]) / 2
+
+    f_avrg = spk.average_profile([f1, f2])
+
+    assert_array_almost_equal(f_avrg.x, x_expected, decimal=16)
+    assert_array_almost_equal(f_avrg.y1, y1_expected, decimal=16)
+    assert_array_almost_equal(f_avrg.y2, y2_expected, decimal=16)
+
+
 if __name__ == "__main__":
     test_pwc()
     test_pwc_add()