diff options
-rw-r--r-- | pyspike/function.py | 106 | ||||
-rw-r--r-- | test/test_function.py | 49 |
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() |