The integral() function in PieceWiseLinFunc.py returns an incorrect (potentially negative) value when the interval argument (start, end) is such that there are no spikes between start and end.
Below is a possible fix that handles small intervals as a separate case.
def integral(self, interval=None):
""" Returns the integral over the given interval.
: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 integral
:rtype: float
"""
def intermediate_value(x0, x1, y0, y1, x):
""" computes the intermediate value of a linear function """
return y0 + (y1-y0)*(x-x0)/(x1-x0)
if interval is None:
# 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')
end_ind = np.searchsorted(self.x, interval[1], side='left')-1
assert start_ind > 0 and end_ind < len(self.x), \
"Invalid averaging interval"
if start_ind <= end_ind:
# first the contribution from between the indices
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
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
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]
))
### if we are here, the interval is completely enclosed between two spikes
### with no spikes in between
else:
# first the contribution from between the indices
start_ind, end_ind = end_ind, start_ind
assert start_ind >= 0 and end_ind < len(self.x), \
"Invalid averaging interval"
assert (end_ind - start_ind) == 1
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
integral -= (interval[0]-self.x[start_ind]) * 0.5 * \
(self.y1[start_ind] +
intermediate_value(self.x[start_ind],
self.x[end_ind],
self.y1[start_ind],
self.y2[start_ind],
interval[0]
))
# correction from last index to end
integral -= (self.x[end_ind] - interval[1]) * 0.5 * \
(self.y2[start_ind] +
intermediate_value(self.x[start_ind], self.x[end_ind],
self.y1[start_ind], self.y2[start_ind],
interval[1]
))
return integral
The integral() function in PieceWiseLinFunc.py returns an incorrect (potentially negative) value when the interval argument (start, end) is such that there are no spikes between start and end. Below is a possible fix that handles small intervals as a separate case.