summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--pyspike/function.py106
-rw-r--r--test/test_function.py49
2 files changed, 150 insertions, 5 deletions
diff --git a/pyspike/function.py b/pyspike/function.py
index ef6b0f1..b705293 100644
--- a/pyspike/function.py
+++ b/pyspike/function.py
@@ -10,6 +10,9 @@ from __future__ import print_function
import numpy as np
+##############################################################
+# PieceWiseConstFunc
+##############################################################
class PieceWiseConstFunc:
""" A class representing a piece-wise constant function. """
@@ -106,6 +109,10 @@ class PieceWiseConstFunc:
self.x = x_new[:index+2]
self.y = y_new[:index+1]
+
+##############################################################
+# PieceWiseLinFunc
+##############################################################
class PieceWiseLinFunc:
""" A class representing a piece-wise linear function. """
@@ -119,9 +126,9 @@ class PieceWiseLinFunc:
- y2: array of length N defining the function values at the right of the
intervals.
"""
- self.x = x
- self.y1 = y1
- self.y2 = y2
+ self.x = np.array(x)
+ self.y1 = np.array(y1)
+ self.y2 = np.array(y2)
def get_plottable_data(self):
""" Returns two arrays containing x- and y-coordinates for immeditate
@@ -136,3 +143,96 @@ class PieceWiseLinFunc:
y_plot[1::2] = self.y2
return x_plot, y_plot
+ def avrg(self):
+ """ Computes the average of the piece-wise linear 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]) * 0.5*(self.y1+self.y2)) / \
+ (self.x[-1]-self.x[0])
+
+ def abs_avrg(self):
+ """ Computes the absolute average of the piece-wise linear 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]) * 0.5 *
+ (np.abs(self.y1)+np.abs(self.y2)))/(self.x[-1]-self.x[0])
+
+ def add(self, f):
+ """ Adds another PieceWiseLin function to this function.
+ Note: only functions defined on the same interval can be summed.
+ Params:
+ - f: PieceWiseLin function to be added.
+ """
+ assert self.x[0] == f.x[0], "The functions have different intervals"
+ assert self.x[-1] == f.x[-1], "The functions have different intervals"
+ x_new = np.empty(len(self.x) + len(f.x))
+ y1_new = np.empty(len(x_new)-1)
+ y2_new = np.empty_like(y1_new)
+ x_new[0] = self.x[0]
+ y1_new[0] = self.y1[0] + f.y1[0]
+ index1 = 0 # index for self
+ index2 = 0 # index for f
+ index = 0 # index for new
+ while (index1+1 < len(self.y1)) and (index2+1 < len(f.y1)):
+ # print(index1+1, self.x[index1+1], self.y[index1+1], x_new[index])
+ if self.x[index1+1] < f.x[index2+1]:
+ # first compute the end value of the previous interval
+ # linear interpolation of the interval
+ y = f.y1[index2] + (f.y2[index2]-f.y1[index2]) * \
+ (self.x[index1+1]-f.x[index2]) / (f.x[index2+1]-f.x[index2])
+ y2_new[index] = self.y2[index1] + y
+ index1 += 1
+ index += 1
+ x_new[index] = self.x[index1]
+ # and the starting value for the next interval
+ y1_new[index] = self.y1[index1] + y
+ elif self.x[index1+1] > f.x[index2+1]:
+ # first compute the end value of the previous interval
+ # linear interpolation of the interval
+ y = self.y1[index1] + (self.y2[index1]-self.y1[index1]) * \
+ (f.x[index2+1]-self.x[index1]) / \
+ (self.x[index1+1]-self.x[index1])
+ y2_new[index] = f.y2[index2] + y
+ index2 += 1
+ index += 1
+ x_new[index] = f.x[index2]
+ # and the starting value for the next interval
+ y1_new[index] = f.y1[index2] + y
+ else: # self.x[index1+1] == f.x[index2+1]:
+ y2_new[index] = self.y2[index1] + f.y2[index2]
+ index1 += 1
+ index2 += 1
+ index += 1
+ x_new[index] = self.x[index1]
+ y1_new[index] = self.y1[index1] + f.y1[index2]
+ # one array reached the end -> copy the contents of the other to the end
+ if index1+1 < len(self.y1):
+ # compute the linear interpolations values
+ y = f.y1[index2] + (f.y2[index2]-f.y1[index2]) * \
+ (self.x[index1+1:-1]-f.x[index2]) / (f.x[index2+1]-f.x[index2])
+ x_new[index+1:index+1+len(self.x)-index1-1] = self.x[index1+1:]
+ y1_new[index+1:index+1+len(self.y1)-index1-1] = self.y1[index1+1:]+y
+ y2_new[index:index+len(self.y2)-index1-1] = self.y2[index1:-1] + y
+ index += len(self.x)-index1-2
+ elif index2+1 < len(f.y1):
+ # compute the linear interpolations values
+ y = self.y1[index1] + (self.y2[index1]-self.y1[index1]) * \
+ (f.x[index2+1:-1]-self.x[index1]) / \
+ (self.x[index1+1]-self.x[index1])
+ x_new[index+1:index+1+len(f.x)-index2-1] = f.x[index2+1:]
+ y1_new[index+1:index+1+len(f.y1)-index2-1] = f.y1[index2+1:] + y
+ y2_new[index:index+len(f.y2)-index2-1] = f.y2[index2:-1] + y
+ index += len(f.x)-index2-2
+ else: # both arrays reached the end simultaneously
+ # only the last x-value missing
+ x_new[index+1] = self.x[-1]
+ # finally, the end value for the last interval
+ y2_new[index] = self.y2[-1]+f.y2[-1]
+ # only use the data that was actually filled
+ self.x = x_new[:index+2]
+ self.y1 = y1_new[:index+1]
+ self.y2 = y2_new[:index+1]
diff --git a/test/test_function.py b/test/test_function.py
index 386b999..014ecac 100644
--- a/test/test_function.py
+++ b/test/test_function.py
@@ -22,8 +22,8 @@ def test_pwc():
xp_expected = [0.0, 1.0, 1.0, 2.0, 2.0, 2.5, 2.5, 4.0]
yp_expected = [1.0, 1.0, -0.5, -0.5, 1.5, 1.5, 0.75, 0.75]
- assert_array_almost_equal(xp, xp_expected)
- assert_array_almost_equal(yp, yp_expected)
+ assert_array_almost_equal(xp, xp_expected, decimal=16)
+ assert_array_almost_equal(yp, yp_expected, decimal=16)
assert_almost_equal(f.avrg(), (1.0-0.5+0.5*1.5+1.5*0.75)/4.0, decimal=16)
assert_almost_equal(f.abs_avrg(), (1.0+0.5+0.5*1.5+1.5*0.75)/4.0,
@@ -49,5 +49,50 @@ def test_pwc():
assert_array_almost_equal(f1.y, 2*f2.y, 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]
+ y2 = [1.5, -0.4, 1.5, 0.25]
+ f = spk.PieceWiseLinFunc(x, y1, y2)
+ xp, yp = f.get_plottable_data()
+
+ xp_expected = [0.0, 1.0, 1.0, 2.0, 2.0, 2.5, 2.5, 4.0]
+ yp_expected = [1.0, 1.5, -0.5, -0.4, 1.5, 1.5, 0.75, 0.25]
+ assert_array_almost_equal(xp, xp_expected, decimal=16)
+ assert_array_almost_equal(yp, yp_expected, decimal=16)
+
+ avrg_expected = (1.25 - 0.45 + 0.75 + 1.5*0.5) / 4.0
+ assert_almost_equal(f.avrg(), avrg_expected, decimal=16)
+
+ abs_avrg_expected = (1.25 + 0.45 + 0.75 + 1.5*0.5) / 4.0
+ assert_almost_equal(f.abs_avrg(), abs_avrg_expected, decimal=16)
+
+ f1 = copy(f)
+ 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)
+ f1.add(f2)
+ x_expected = [0.0, 0.75, 1.0, 2.0, 2.5, 2.7, 4.0]
+ y1_expected = [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]
+ y2_expected = [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]
+ assert_array_almost_equal(f1.x, x_expected, decimal=16)
+ assert_array_almost_equal(f1.y1, y1_expected, decimal=16)
+ assert_array_almost_equal(f1.y2, y2_expected, decimal=16)
+
+ f2.add(f)
+ assert_array_almost_equal(f2.x, x_expected, decimal=16)
+ assert_array_almost_equal(f2.y1, y1_expected, decimal=16)
+ assert_array_almost_equal(f2.y2, y2_expected, decimal=16)
+
+ f1.add(f2)
+ # same x, but y doubled
+ assert_array_almost_equal(f1.x, f2.x, decimal=16)
+ assert_array_almost_equal(f1.y1, 2*f2.y1, decimal=16)
+ assert_array_almost_equal(f1.y2, 2*f2.y2, decimal=16)
+
+
if __name__ == "__main__":
test_pwc()