diff options
author | Mario Mulansky <mario.mulansky@gmx.net> | 2014-10-23 15:33:06 +0200 |
---|---|---|
committer | Mario Mulansky <mario.mulansky@gmx.net> | 2014-10-23 15:33:06 +0200 |
commit | 3a36f81d52137435910a4b7c656a478ae5b38ece (patch) | |
tree | f44c3f31e8fdf598e1b406c4ada24bcb2280c1cc | |
parent | ac8e9ca85a889ee2d9fe984d19976917ffdfea46 (diff) |
support for multiple averaging intervals
-rw-r--r-- | examples/averages.py | 26 | ||||
-rw-r--r-- | pyspike/function.py | 144 | ||||
-rw-r--r-- | test/test_function.py | 10 |
3 files changed, 126 insertions, 54 deletions
diff --git a/examples/averages.py b/examples/averages.py index a06871a..6413ab2 100644 --- a/examples/averages.py +++ b/examples/averages.py @@ -20,23 +20,25 @@ print("ISI-distance: %.8f" % f.avrg()) isi1 = f.avrg(interval=(0, 1000)) isi2 = f.avrg(interval=(1000, 2000)) -isi3 = f.avrg(interval=(2000, 3000)) -isi4 = f.avrg(interval=(3000, 4000)) +isi3 = f.avrg(interval=[(0, 1000), (2000, 3000)]) +isi4 = f.avrg(interval=[(1000, 2000), (3000, 4000)]) -print("ISI-distance (0-1000): %.8f" % isi1) -print("ISI-distance (1000-2000): %.8f" % isi2) -print("ISI-distance (2000-3000): %.8f" % isi3) -print("ISI-distance (3000-4000): %.8f" % isi4) +print("ISI-distance (0-1000): %.8f" % isi1) +print("ISI-distance (1000-2000): %.8f" % isi2) +print("ISI-distance (0-1000) and (2000-3000): %.8f" % isi3) +print("ISI-distance (1000-2000) and (3000-4000): %.8f" % isi4) print() +f = spk.spike_profile(spike_trains[0], spike_trains[1]) + print("SPIKE-distance: %.8f" % f.avrg()) spike1 = f.avrg(interval=(0, 1000)) spike2 = f.avrg(interval=(1000, 2000)) -spike3 = f.avrg(interval=(2000, 3000)) -spike4 = f.avrg(interval=(3000, 4000)) +spike3 = f.avrg(interval=[(0, 1000), (2000, 3000)]) +spike4 = f.avrg(interval=[(1000, 2000), (3000, 4000)]) -print("SPIKE-distance (0-1000): %.8f" % spike1) -print("SPIKE-distance (1000-2000): %.8f" % spike2) -print("SPIKE-distance (2000-3000): %.8f" % spike3) -print("SPIKE-distance (3000-4000): %.8f" % spike4) +print("SPIKE-distance (0-1000): %.8f" % spike1) +print("SPIKE-distance (1000-2000): %.8f" % spike2) +print("SPIKE-distance (0-1000) and (2000-3000): %.8f" % spike3) +print("SPIKE-distance (1000-2000) and (3000-4000): %.8f" % spike4) diff --git a/pyspike/function.py b/pyspike/function.py index d53042f..0c0a391 100644 --- a/pyspike/function.py +++ b/pyspike/function.py @@ -11,6 +11,7 @@ Distributed under the BSD License from __future__ import print_function import numpy as np +import collections ############################################################## @@ -74,21 +75,18 @@ class PieceWiseConstFunc(object): return x_plot, y_plot - def avrg(self, interval=None): - """ Computes the average of the piece-wise const function: - :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval. + def integral(self, interval=None): + """ Returns the integral over the given interval. - :param interval: averaging interval given as a pair of floats, if None - the average over the whole function is computed. + :param interval: integration interval given as a pair of floats, if + None the integral over the whole function is computed. :type interval: Pair of floats or None. - :returns: the average a. - :rtype: double - + :returns: the integral + :rtype: float """ if interval is None: - # no interval given, average over the whole spike train - a = np.sum((self.x[1:]-self.x[:-1]) * self.y) / \ - (self.x[-1]-self.x[0]) + # no interval given, integrate over the whole spike train + a = np.sum((self.x[1:]-self.x[:-1]) * self.y) else: # find the indices corresponding to the interval start_ind = np.searchsorted(self.x, interval[0], side='right') @@ -103,7 +101,39 @@ class PieceWiseConstFunc(object): a += (self.x[start_ind]-interval[0]) * self.y[start_ind-1] # correction from last index to end a += (interval[1]-self.x[end_ind]) * self.y[end_ind] - a /= (interval[1]-interval[0]) + return a + + def avrg(self, interval=None): + """ Computes the average of the piece-wise const function: + :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval. + + :param interval: averaging interval given as a pair of floats, a + sequence of pairs for averaging multiple intervals, or + None, if None the average over the whole function is + computed. + :type interval: Pair, sequence of pairs, or None. + :returns: the average a. + :rtype: float + """ + if interval is None: + # no interval given, average over the whole spike train + return self.integral() / (self.x[-1]-self.x[0]) + + # check if interval is as sequence + assert isinstance(interval, collections.Sequence), \ + "Invalid value for `interval`. None, Sequence or Tuple expected." + # check if interval is a sequence of intervals + if not isinstance(interval[0], collections.Sequence): + # just one interval + a = self.integral(interval) / (interval[1]-interval[0]) + else: + # several intervals + a = 0.0 + int_length = 0.0 + for ival in interval: + a += self.integral(ival) + int_length += ival[1] - ival[0] + a /= int_length return a def add(self, f): @@ -203,16 +233,14 @@ class PieceWiseLinFunc: y_plot[1::2] = self.y2 return x_plot, y_plot - def avrg(self, interval=None): - """ Computes the average of the piece-wise linear function: - :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval. + def integral(self, interval=None): + """ Returns the integral over the given interval. - :param interval: averaging interval given as a pair of floats, if None - the average over the whole function is computed. + :param interval: integration interval given as a pair of floats, if + None the integral over the whole function is computed. :type interval: Pair of floats or None. - :returns: the average a. - :rtype: double - + :returns: the integral + :rtype: float """ def intermediate_value(x0, x1, y0, y1, x): @@ -220,9 +248,8 @@ class PieceWiseLinFunc: return y0 + (y1-y0)*(x-x0)/(x1-x0) if interval is None: - # no interval given, average over the whole spike train - a = np.sum((self.x[1:]-self.x[:-1]) * 0.5*(self.y1+self.y2)) / \ - (self.x[-1]-self.x[0]) + # no interval given, integrate over the whole spike train + integral = np.sum((self.x[1:]-self.x[:-1]) * 0.5*(self.y1+self.y2)) else: # find the indices corresponding to the interval start_ind = np.searchsorted(self.x, interval[0], side='right') @@ -230,26 +257,61 @@ class PieceWiseLinFunc: assert start_ind > 0 and end_ind < len(self.x), \ "Invalid averaging interval" # first the contribution from between the indices - a = np.sum((self.x[start_ind+1:end_ind+1] - - self.x[start_ind:end_ind]) * - 0.5*(self.y1[start_ind:end_ind] + - self.y2[start_ind:end_ind])) + integral = np.sum((self.x[start_ind+1:end_ind+1] - + self.x[start_ind:end_ind]) * + 0.5*(self.y1[start_ind:end_ind] + + self.y2[start_ind:end_ind])) # correction from start to first index - a += (self.x[start_ind]-interval[0]) * 0.5 * \ - (self.y2[start_ind-1] + - intermediate_value(self.x[start_ind-1], self.x[start_ind], - self.y1[start_ind-1], - self.y2[start_ind-1], - interval[0] - )) + integral += (self.x[start_ind]-interval[0]) * 0.5 * \ + (self.y2[start_ind-1] + + intermediate_value(self.x[start_ind-1], + self.x[start_ind], + self.y1[start_ind-1], + self.y2[start_ind-1], + interval[0] + )) # correction from last index to end - a += (interval[1]-self.x[end_ind]) * 0.5 * \ - (self.y1[end_ind] + - intermediate_value(self.x[end_ind], self.x[end_ind+1], - self.y1[end_ind], self.y2[end_ind], - interval[1] - )) - a /= (interval[1]-interval[0]) + integral += (interval[1]-self.x[end_ind]) * 0.5 * \ + (self.y1[end_ind] + + intermediate_value(self.x[end_ind], self.x[end_ind+1], + self.y1[end_ind], self.y2[end_ind], + interval[1] + )) + return integral + + def avrg(self, interval=None): + """ Computes the average of the piece-wise linear function: + :math:`a = 1/T int_0^T f(x) dx` where T is the length of the interval. + + :param interval: averaging interval given as a pair of floats, a + sequence of pairs for averaging multiple intervals, or + None, if None the average over the whole function is + computed. + :type interval: Pair, sequence of pairs, or None. + :returns: the average a. + :rtype: float + + """ + + if interval is None: + # no interval given, average over the whole spike train + return self.integral() / (self.x[-1]-self.x[0]) + + # check if interval is as sequence + assert isinstance(interval, collections.Sequence), \ + "Invalid value for `interval`. None, Sequence or Tuple expected." + # check if interval is a sequence of intervals + if not isinstance(interval[0], collections.Sequence): + # just one interval + a = self.integral(interval) / (interval[1]-interval[0]) + else: + # several intervals + a = 0.0 + int_length = 0.0 + for ival in interval: + a += self.integral(ival) + int_length += ival[1] - ival[0] + a /= int_length return a def add(self, f): diff --git a/test/test_function.py b/test/test_function.py index ed0b2ed..ba87db8 100644 --- a/test/test_function.py +++ b/test/test_function.py @@ -39,6 +39,10 @@ def test_pwc(): a = f.avrg([1.0, 4.0]) assert_almost_equal(a, (-0.5*1.0+0.5*1.5+1.5*0.75)/3.0, decimal=16) + # averaging over multiple intervals + a = f.avrg([(0.5, 1.5), (1.5, 3.5)]) + assert_almost_equal(a, (0.5-0.5+0.5*1.5+1.0*0.75)/3.0, decimal=16) + def test_pwc_add(): # some random data @@ -116,12 +120,16 @@ def test_pwl(): a = f.avrg([1.5, 3.5]) assert_almost_equal(a, (-0.425*0.5 + 0.75 + (0.75+0.75-0.5/1.5)/2) / 2.0, decimal=16) - a = f.avrg([1.0, 3.5]) + a = f.avrg((1.0, 3.5)) assert_almost_equal(a, (-0.45 + 0.75 + (0.75+0.75-0.5/1.5)/2) / 2.5, decimal=16) a = f.avrg([1.0, 4.0]) assert_almost_equal(a, (-0.45 + 0.75 + 1.5*0.5) / 3.0, decimal=16) + # averaging over multiple intervals + a = f.avrg([(0.5, 1.5), (1.5, 2.5)]) + assert_almost_equal(a, (1.375*0.5 - 0.45 + 0.75)/2.0, decimal=16) + def test_pwl_add(): x = [0.0, 1.0, 2.0, 2.5, 4.0] |