summaryrefslogtreecommitdiff
path: root/pyspike
diff options
context:
space:
mode:
Diffstat (limited to 'pyspike')
-rw-r--r--pyspike/DiscreteFunc.py21
-rw-r--r--pyspike/PieceWiseConstFunc.py46
-rw-r--r--pyspike/PieceWiseLinFunc.py58
-rw-r--r--pyspike/SpikeTrain.py10
-rw-r--r--pyspike/cython/cython_add.pyx20
-rw-r--r--pyspike/cython/cython_distances.pyx398
-rw-r--r--pyspike/cython/cython_profiles.pyx (renamed from pyspike/cython/cython_distance.pyx)51
-rw-r--r--pyspike/cython/python_backend.py14
-rw-r--r--pyspike/generic.py77
-rw-r--r--pyspike/isi_distance.py39
-rw-r--r--pyspike/spike_distance.py44
-rw-r--r--pyspike/spike_sync.py79
12 files changed, 765 insertions, 92 deletions
diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py
index 33b7a81..17153ee 100644
--- a/pyspike/DiscreteFunc.py
+++ b/pyspike/DiscreteFunc.py
@@ -125,17 +125,21 @@ class DiscreteFunc(object):
def integral(self, interval=None):
""" Returns the integral over the given interval. For the discrete
- function, this amounts to the sum over all values divided by the total
- multiplicity.
+ function, this amounts to two values: the sum over all values and the
+ sum over all multiplicities.
:param interval: integration interval given as a pair of floats, or a
sequence of pairs in case of multiple intervals, if
None the integral over the whole function is computed.
:type interval: Pair, sequence of pairs, or None.
- :returns: the integral
- :rtype: float
+ :returns: the summed values and the summed multiplicity
+ :rtype: pair of float
"""
+ if len(self.y) <= 2:
+ # no actual values in the profile, return spike sync of 1
+ return 1.0, 1.0
+
def get_indices(ival):
""" Retuns the indeces surrounding the given interval"""
start_ind = np.searchsorted(self.x, ival[0], side='right')
@@ -147,7 +151,7 @@ class DiscreteFunc(object):
if interval is None:
# no interval given, integrate over the whole spike train
# don't count the first value, which is zero by definition
- return 1.0 * np.sum(self.y[1:-1]) / np.sum(self.mp[1:-1])
+ return (1.0 * np.sum(self.y[1:-1]), np.sum(self.mp[1:-1]))
# check if interval is as sequence
assert isinstance(interval, collections.Sequence), \
@@ -156,7 +160,7 @@ class DiscreteFunc(object):
if not isinstance(interval[0], collections.Sequence):
# find the indices corresponding to the interval
start_ind, end_ind = get_indices(interval)
- return (np.sum(self.y[start_ind:end_ind]) /
+ return (np.sum(self.y[start_ind:end_ind]),
np.sum(self.mp[start_ind:end_ind]))
else:
value = 0.0
@@ -166,7 +170,7 @@ class DiscreteFunc(object):
start_ind, end_ind = get_indices(ival)
value += np.sum(self.y[start_ind:end_ind])
multiplicity += np.sum(self.mp[start_ind:end_ind])
- return value/multiplicity
+ return (value, multiplicity)
def avrg(self, interval=None):
""" Computes the average of the interval sequence:
@@ -180,7 +184,8 @@ class DiscreteFunc(object):
:returns: the average a.
:rtype: float
"""
- return self.integral(interval)
+ val, mp = self.integral(interval)
+ return val/mp
def add(self, f):
""" Adds another `DiscreteFunc` function to this function.
diff --git a/pyspike/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py
index 41998ef..2705443 100644
--- a/pyspike/PieceWiseConstFunc.py
+++ b/pyspike/PieceWiseConstFunc.py
@@ -26,6 +26,52 @@ class PieceWiseConstFunc(object):
self.x = np.array(x)
self.y = np.array(y)
+ def __call__(self, t):
+ """ Returns the function value for the given time t. If t is a list of
+ times, the corresponding list of values is returned.
+
+ :param: time t, or list of times
+ :returns: function value(s) at that time(s).
+ """
+ assert np.all(t >= self.x[0]) and np.all(t <= self.x[-1]), \
+ "Invalid time: " + str(t)
+
+ ind = np.searchsorted(self.x, t, side='right')
+
+ if isinstance(t, collections.Sequence):
+ # t is a sequence of values
+ # correct the cases t == x[0], t == x[-1]
+ ind[ind == 0] = 1
+ ind[ind == len(self.x)] = len(self.x)-1
+ value = self.y[ind-1]
+ # correct the values at exact spike times: there the value should
+ # be the at half of the step
+ # obtain the 'left' side indices for t
+ ind_l = np.searchsorted(self.x, t, side='left')
+ # if left and right side indices differ, the time t has to appear
+ # in self.x
+ ind_at_spike = np.logical_and(np.logical_and(ind != ind_l,
+ ind > 1),
+ ind < len(self.x))
+ # get the corresponding indices for the resulting value array
+ val_ind = np.arange(len(ind))[ind_at_spike]
+ # and for the arrays self.x, y1, y2
+ xy_ind = ind[ind_at_spike]
+ # use the middle of the left and right ISI value
+ value[val_ind] = 0.5 * (self.y[xy_ind-1] + self.y[xy_ind-2])
+ return value
+ else: # t is a single value
+ # specific check for interval edges
+ if t == self.x[0]:
+ return self.y[0]
+ if t == self.x[-1]:
+ return self.y[-1]
+ # check if we are on any other exact spike time
+ if sum(self.x == t) > 0:
+ # use the middle of the left and right ISI value
+ return 0.5 * (self.y[ind-1] + self.y[ind-2])
+ return self.y[ind-1]
+
def copy(self):
""" Returns a copy of itself
diff --git a/pyspike/PieceWiseLinFunc.py b/pyspike/PieceWiseLinFunc.py
index f2442be..c0dd475 100644
--- a/pyspike/PieceWiseLinFunc.py
+++ b/pyspike/PieceWiseLinFunc.py
@@ -29,6 +29,64 @@ class PieceWiseLinFunc:
self.y1 = np.array(y1)
self.y2 = np.array(y2)
+ def __call__(self, t):
+ """ Returns the function value for the given time t. If t is a list of
+ times, the corresponding list of values is returned.
+
+ :param: time t, or list of times
+ :returns: function value(s) at that time(s).
+ """
+ def intermediate_value(x0, x1, y0, y1, x):
+ """ computes the intermediate value of a linear function """
+ return y0 + (y1-y0)*(x-x0)/(x1-x0)
+
+ assert np.all(t >= self.x[0]) and np.all(t <= self.x[-1]), \
+ "Invalid time: " + str(t)
+
+ ind = np.searchsorted(self.x, t, side='right')
+
+ if isinstance(t, collections.Sequence):
+ # t is a sequence of values
+ # correct the cases t == x[0], t == x[-1]
+ ind[ind == 0] = 1
+ ind[ind == len(self.x)] = len(self.x)-1
+ value = intermediate_value(self.x[ind-1],
+ self.x[ind],
+ self.y1[ind-1],
+ self.y2[ind-1],
+ t)
+ # correct the values at exact spike times: there the value should
+ # be the at half of the step
+ # obtain the 'left' side indices for t
+ ind_l = np.searchsorted(self.x, t, side='left')
+ # if left and right side indices differ, the time t has to appear
+ # in self.x
+ ind_at_spike = np.logical_and(np.logical_and(ind != ind_l,
+ ind > 1),
+ ind < len(self.x))
+ # get the corresponding indices for the resulting value array
+ val_ind = np.arange(len(ind))[ind_at_spike]
+ # and for the values in self.x, y1, y2
+ xy_ind = ind[ind_at_spike]
+ # the values are defined as the average of the left and right limit
+ value[val_ind] = 0.5 * (self.y1[xy_ind-1] + self.y2[xy_ind-2])
+ return value
+ else: # t is a single value
+ # specific check for interval edges
+ if t == self.x[0]:
+ return self.y1[0]
+ if t == self.x[-1]:
+ return self.y2[-1]
+ # check if we are on any other exact spike time
+ if sum(self.x == t) > 0:
+ # use the middle of the left and right Spike value
+ return 0.5 * (self.y1[ind-1] + self.y2[ind-2])
+ return intermediate_value(self.x[ind-1],
+ self.x[ind],
+ self.y1[ind-1],
+ self.y2[ind-1],
+ t)
+
def copy(self):
""" Returns a copy of itself
diff --git a/pyspike/SpikeTrain.py b/pyspike/SpikeTrain.py
index a02b7ab..9127b60 100644
--- a/pyspike/SpikeTrain.py
+++ b/pyspike/SpikeTrain.py
@@ -48,3 +48,13 @@ class SpikeTrain(object):
"""
return SpikeTrain(self.spikes.copy(), [self.t_start, self.t_end])
+
+ def get_spikes_non_empty(self):
+ """Returns the spikes of this spike train with auxiliary spikes in case
+ of empty spike trains.
+ """
+ if len(self.spikes) < 2:
+ return np.unique(np.insert([self.t_start, self.t_end], 1,
+ self.spikes))
+ else:
+ return self.spikes
diff --git a/pyspike/cython/cython_add.pyx b/pyspike/cython/cython_add.pyx
index ac64005..8da1e53 100644
--- a/pyspike/cython/cython_add.pyx
+++ b/pyspike/cython/cython_add.pyx
@@ -83,13 +83,9 @@ def add_piece_wise_const_cython(double[:] x1, double[:] y1,
else: # both arrays reached the end simultaneously
# only the last x-value missing
x_new[index+1] = x1[N1-1]
- # the last value is again the end of the interval
- # x_new[index+1] = x1[-1]
- # only use the data that was actually filled
- x1 = x_new[:index+2]
- y1 = y_new[:index+1]
# end nogil
- return np.array(x_new[:index+2]), np.array(y_new[:index+1])
+ # return np.asarray(x_new[:index+2]), np.asarray(y_new[:index+1])
+ return np.asarray(x_new[:index+2]), np.asarray(y_new[:index+1])
############################################################
@@ -169,9 +165,9 @@ def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12,
y2_new[index] = y12[N1-2]+y22[N2-2]
# only use the data that was actually filled
# end nogil
- return (np.array(x_new[:index+2]),
- np.array(y1_new[:index+1]),
- np.array(y2_new[:index+1]))
+ return (np.asarray(x_new[:index+2]),
+ np.asarray(y1_new[:index+1]),
+ np.asarray(y2_new[:index+1]))
############################################################
@@ -230,6 +226,6 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1,
# the last value is again the end of the interval
# only use the data that was actually filled
- return (np.array(x_new[:index+1]),
- np.array(y_new[:index+1]),
- np.array(mp_new[:index+1]))
+ return (np.asarray(x_new[:index+1]),
+ np.asarray(y_new[:index+1]),
+ np.asarray(mp_new[:index+1]))
diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx
new file mode 100644
index 0000000..c4f2349
--- /dev/null
+++ b/pyspike/cython/cython_distances.pyx
@@ -0,0 +1,398 @@
+#cython: boundscheck=False
+#cython: wraparound=False
+#cython: cdivision=True
+
+"""
+cython_distances.pyx
+
+cython implementation of the isi-, spike- and spike-sync distances
+
+Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects
+improves the performance of spike_distance by a factor of 10!
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+
+"""
+To test whether things can be optimized: remove all yellow stuff
+in the html output::
+
+ cython -a cython_distances.pyx
+
+which gives::
+
+ cython_distances.html
+
+"""
+
+import numpy as np
+cimport numpy as np
+
+from libc.math cimport fabs
+from libc.math cimport fmax
+from libc.math cimport fmin
+
+DTYPE = np.float
+ctypedef np.float_t DTYPE_t
+
+
+############################################################
+# isi_distance_cython
+############################################################
+def isi_distance_cython(double[:] s1, double[:] s2,
+ double t_start, double t_end):
+
+ cdef double isi_value
+ cdef int index1, index2, index
+ cdef int N1, N2
+ cdef double nu1, nu2
+ cdef double last_t, curr_t, curr_isi
+ isi_value = 0.0
+ N1 = len(s1)
+ N2 = len(s2)
+
+ # first interspike interval - check if a spike exists at the start time
+ if s1[0] > t_start:
+ # edge correction
+ nu1 = fmax(s1[0]-t_start, s1[1]-s1[0])
+ index1 = -1
+ else:
+ nu1 = s1[1]-s1[0]
+ index1 = 0
+
+ if s2[0] > t_start:
+ # edge correction
+ nu2 = fmax(s2[0]-t_start, s2[1]-s2[0])
+ index2 = -1
+ else:
+ nu2 = s2[1]-s2[0]
+ index2 = 0
+
+ last_t = t_start
+ curr_isi = fabs(nu1-nu2)/fmax(nu1, nu2)
+ index = 1
+
+ with nogil: # release the interpreter to allow multithreading
+ while index1+index2 < N1+N2-2:
+ # check which spike is next, only if there are spikes left in 1
+ # next spike in 1 is earlier, or there are no spikes left in 2
+ if (index1 < N1-1) and ((index2 == N2-1) or
+ (s1[index1+1] < s2[index2+1])):
+ index1 += 1
+ curr_t = s1[index1]
+ if index1 < N1-1:
+ nu1 = s1[index1+1]-s1[index1]
+ else:
+ # edge correction
+ nu1 = fmax(t_end-s1[index1], nu1)
+ elif (index2 < N2-1) and ((index1 == N1-1) or
+ (s1[index1+1] > s2[index2+1])):
+ index2 += 1
+ curr_t = s2[index2]
+ if index2 < N2-1:
+ nu2 = s2[index2+1]-s2[index2]
+ else:
+ # edge correction
+ nu2 = fmax(t_end-s2[index2], nu2)
+ else: # s1[index1+1] == s2[index2+1]
+ index1 += 1
+ index2 += 1
+ curr_t = s1[index1]
+ if index1 < N1-1:
+ nu1 = s1[index1+1]-s1[index1]
+ else:
+ # edge correction
+ nu1 = fmax(t_end-s1[index1], nu1)
+ if index2 < N2-1:
+ nu2 = s2[index2+1]-s2[index2]
+ else:
+ # edge correction
+ nu2 = fmax(t_end-s2[index2], nu2)
+ # compute the corresponding isi-distance
+ isi_value += curr_isi * (curr_t - last_t)
+ curr_isi = fabs(nu1 - nu2) / fmax(nu1, nu2)
+ last_t = curr_t
+ index += 1
+
+ isi_value += curr_isi * (t_end - last_t)
+ # end nogil
+
+ return isi_value / (t_end-t_start)
+
+
+############################################################
+# get_min_dist_cython
+############################################################
+cdef inline double get_min_dist_cython(double spike_time,
+ double[:] spike_train,
+ # use memory view to ensure inlining
+ # np.ndarray[DTYPE_t,ndim=1] spike_train,
+ int N,
+ int start_index,
+ double t_start, double t_end) nogil:
+ """ Returns the minimal distance |spike_time - spike_train[i]|
+ with i>=start_index.
+ """
+ cdef double d, d_temp
+ # start with the distance to the start time
+ d = fabs(spike_time - t_start)
+ if start_index < 0:
+ start_index = 0
+ while start_index < N:
+ d_temp = fabs(spike_time - spike_train[start_index])
+ if d_temp > d:
+ return d
+ else:
+ d = d_temp
+ start_index += 1
+
+ # finally, check the distance to end time
+ d_temp = fabs(t_end - spike_time)
+ if d_temp > d:
+ return d
+ else:
+ return d_temp
+
+
+############################################################
+# isi_avrg_cython
+############################################################
+cdef inline double isi_avrg_cython(double isi1, double isi2) nogil:
+ return 0.5*(isi1+isi2)*(isi1+isi2)
+ # alternative definition to obtain <S> ~ 0.5 for Poisson spikes
+ # return 0.5*(isi1*isi1+isi2*isi2)
+
+
+############################################################
+# spike_distance_cython
+############################################################
+def spike_distance_cython(double[:] t1, double[:] t2,
+ double t_start, double t_end):
+
+ cdef int N1, N2, index1, index2, index
+ cdef double t_p1, t_f1, t_p2, t_f2, dt_p1, dt_p2, dt_f1, dt_f2
+ cdef double isi1, isi2, s1, s2
+ cdef double y_start, y_end, t_last, t_current, spike_value
+
+ spike_value = 0.0
+
+ N1 = len(t1)
+ N2 = len(t2)
+
+ with nogil: # release the interpreter to allow multithreading
+ t_last = t_start
+ t_p1 = t_start
+ t_p2 = t_start
+ if t1[0] > t_start:
+ # dt_p1 = t2[0]-t_start
+ t_f1 = t1[0]
+ dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end)
+ isi1 = fmax(t_f1-t_start, t1[1]-t1[0])
+ dt_p1 = dt_f1
+ s1 = dt_p1*(t_f1-t_start)/isi1
+ index1 = -1
+ else:
+ t_f1 = t1[1]
+ dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, t_end)
+ dt_p1 = 0.0
+ isi1 = t1[1]-t1[0]
+ s1 = dt_p1
+ index1 = 0
+ if t2[0] > t_start:
+ # dt_p1 = t2[0]-t_start
+ t_f2 = t2[0]
+ dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end)
+ dt_p2 = dt_f2
+ isi2 = fmax(t_f2-t_start, t2[1]-t2[0])
+ s2 = dt_p2*(t_f2-t_start)/isi2
+ index2 = -1
+ else:
+ t_f2 = t2[1]
+ dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, t_end)
+ dt_p2 = 0.0
+ isi2 = t2[1]-t2[0]
+ s2 = dt_p2
+ index2 = 0
+
+ y_start = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+ index = 1
+
+ while index1+index2 < N1+N2-2:
+ # print(index, index1, index2)
+ if (index1 < N1-1) and (t_f1 < t_f2 or index2 == N2-1):
+ index1 += 1
+ # first calculate the previous interval end value
+ s1 = dt_f1*(t_f1-t_p1) / isi1
+ # the previous time now was the following time before:
+ dt_p1 = dt_f1
+ t_p1 = t_f1 # t_p1 contains the current time point
+ # get the next time
+ if index1 < N1-1:
+ t_f1 = t1[index1+1]
+ else:
+ t_f1 = t_end
+ t_curr = t_p1
+ s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
+ y_end = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+
+ spike_value += 0.5*(y_start + y_end) * (t_curr - t_last)
+
+ # now the next interval start value
+ if index1 < N1-1:
+ dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
+ t_start, t_end)
+ isi1 = t_f1-t_p1
+ s1 = dt_p1
+ else:
+ dt_f1 = dt_p1
+ isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ # s1 needs adjustment due to change of isi1
+ s1 = dt_p1*(t_end-t1[N1-1])/isi1
+ # s2 is the same as above, thus we can compute y2 immediately
+ y_start = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+ elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1):
+ index2 += 1
+ # first calculate the previous interval end value
+ s2 = dt_f2*(t_f2-t_p2) / isi2
+ # the previous time now was the following time before:
+ dt_p2 = dt_f2
+ t_p2 = t_f2 # t_p2 contains the current time point
+ # get the next time
+ if index2 < N2-1:
+ t_f2 = t2[index2+1]
+ else:
+ t_f2 = t_end
+ t_curr = t_p2
+ s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
+ y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+
+ spike_value += 0.5*(y_start + y_end) * (t_curr - t_last)
+
+ # now the next interval start value
+ if index2 < N2-1:
+ dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
+ t_start, t_end)
+ isi2 = t_f2-t_p2
+ s2 = dt_p2
+ else:
+ dt_f2 = dt_p2
+ isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ # s2 needs adjustment due to change of isi2
+ s2 = dt_p2*(t_end-t2[N2-1])/isi2
+ # s1 is the same as above, thus we can compute y2 immediately
+ y_start = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
+ else: # t_f1 == t_f2 - generate only one event
+ index1 += 1
+ index2 += 1
+ t_p1 = t_f1
+ t_p2 = t_f2
+ dt_p1 = 0.0
+ dt_p2 = 0.0
+ t_curr = t_f1
+ y_end = 0.0
+ spike_value += 0.5*(y_start + y_end) * (t_curr - t_last)
+ y_start = 0.0
+ if index1 < N1-1:
+ t_f1 = t1[index1+1]
+ dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
+ t_start, t_end)
+ isi1 = t_f1 - t_p1
+ else:
+ t_f1 = t_end
+ dt_f1 = dt_p1
+ isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ if index2 < N2-1:
+ t_f2 = t2[index2+1]
+ dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
+ t_start, t_end)
+ isi2 = t_f2 - t_p2
+ else:
+ t_f2 = t_end
+ dt_f2 = dt_p2
+ isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ index += 1
+ t_last = t_curr
+ # isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ # isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ s1 = dt_f1*(t_end-t1[N1-1])/isi1
+ s2 = dt_f2*(t_end-t2[N2-1])/isi2
+ y_end = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
+ spike_value += 0.5*(y_start + y_end) * (t_end - t_last)
+ # end nogil
+
+ # use only the data added above
+ # could be less than original length due to equal spike times
+ return spike_value / (t_end-t_start)
+
+
+
+############################################################
+# get_tau
+############################################################
+cdef inline double get_tau(double[:] spikes1, double[:] spikes2,
+ int i, int j, double interval, double max_tau):
+ cdef double m = interval # use interval length as initial tau
+ cdef int N1 = spikes1.shape[0]-1 # len(spikes1)-1
+ cdef int N2 = spikes2.shape[0]-1 # len(spikes2)-1
+ if i < N1 and i > -1:
+ m = fmin(m, spikes1[i+1]-spikes1[i])
+ if j < N2 and j > -1:
+ m = fmin(m, spikes2[j+1]-spikes2[j])
+ if i > 0:
+ m = fmin(m, spikes1[i]-spikes1[i-1])
+ if j > 0:
+ m = fmin(m, spikes2[j]-spikes2[j-1])
+ m *= 0.5
+ if max_tau > 0.0:
+ m = fmin(m, max_tau)
+ return m
+
+
+############################################################
+# coincidence_value_cython
+############################################################
+def coincidence_value_cython(double[:] spikes1, double[:] spikes2,
+ double t_start, double t_end, double max_tau):
+
+ cdef int N1 = len(spikes1)
+ cdef int N2 = len(spikes2)
+ cdef int i = -1
+ cdef int j = -1
+ cdef double coinc = 0.0
+ cdef double mp = 0.0
+ cdef double interval = t_end - t_start
+ cdef double tau
+ while i + j < N1 + N2 - 2:
+ if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]):
+ i += 1
+ mp += 1
+ tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
+ if j > -1 and spikes1[i]-spikes2[j] < tau:
+ # coincidence between the current spike and the previous spike
+ # both get marked with 1
+ coinc += 2
+ elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]):
+ j += 1
+ mp += 1
+ tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
+ if i > -1 and spikes2[j]-spikes1[i] < tau:
+ # coincidence between the current spike and the previous spike
+ # both get marked with 1
+ coinc += 2
+ else: # spikes1[i+1] = spikes2[j+1]
+ # advance in both spike trains
+ j += 1
+ i += 1
+ # add only one event, but with coincidence 2 and multiplicity 2
+ mp += 2
+ coinc += 2
+
+ if coinc == 0 and mp == 0:
+ # empty spike trains -> spike sync = 1 by definition
+ coinc = 1
+ mp = 1
+
+ return coinc, mp
diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_profiles.pyx
index 6ee0181..f9893eb 100644
--- a/pyspike/cython/cython_distance.pyx
+++ b/pyspike/cython/cython_profiles.pyx
@@ -3,14 +3,14 @@
#cython: cdivision=True
"""
-cython_distances.pyx
+cython_profiles.pyx
-cython implementation of the isi- and spike-distance
+cython implementation of the isi-, spike- and spike-sync profiles
Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects
improves the performance of spike_distance by a factor of 10!
-Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net>
+Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
Distributed under the BSD License
@@ -20,11 +20,11 @@ Distributed under the BSD License
To test whether things can be optimized: remove all yellow stuff
in the html output::
- cython -a cython_distance.pyx
+ cython -a cython_profiles.pyx
which gives::
- cython_distance.html
+ cython_profiles.html
"""
@@ -40,10 +40,10 @@ ctypedef np.float_t DTYPE_t
############################################################
-# isi_distance_cython
+# isi_profile_cython
############################################################
-def isi_distance_cython(double[:] s1, double[:] s2,
- double t_start, double t_end):
+def isi_profile_cython(double[:] s1, double[:] s2,
+ double t_start, double t_end):
cdef double[:] spike_events
cdef double[:] isi_values
@@ -173,10 +173,10 @@ cdef inline double isi_avrg_cython(double isi1, double isi2) nogil:
############################################################
-# spike_distance_cython
+# spike_profile_cython
############################################################
-def spike_distance_cython(double[:] t1, double[:] t2,
- double t_start, double t_end):
+def spike_profile_cython(double[:] t1, double[:] t2,
+ double t_start, double t_end):
cdef double[:] spike_events
cdef double[:] y_starts
@@ -342,11 +342,11 @@ def spike_distance_cython(double[:] t1, double[:] t2,
############################################################
-# coincidence_python
+# get_tau
############################################################
cdef inline double get_tau(double[:] spikes1, double[:] spikes2,
- int i, int j, double max_tau):
- cdef double m = 1E100 # some huge number
+ int i, int j, double interval, double max_tau):
+ cdef double m = interval # use interval as initial tau
cdef int N1 = spikes1.shape[0]-1 # len(spikes1)-1
cdef int N2 = spikes2.shape[0]-1 # len(spikes2)-1
if i < N1 and i > -1:
@@ -364,10 +364,10 @@ cdef inline double get_tau(double[:] spikes1, double[:] spikes2,
############################################################
-# coincidence_cython
+# coincidence_profile_cython
############################################################
-def coincidence_cython(double[:] spikes1, double[:] spikes2,
- double t_start, double t_end, double max_tau):
+def coincidence_profile_cython(double[:] spikes1, double[:] spikes2,
+ double t_start, double t_end, double max_tau):
cdef int N1 = len(spikes1)
cdef int N2 = len(spikes2)
@@ -377,12 +377,13 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2,
cdef double[:] st = np.zeros(N1 + N2 + 2) # spike times
cdef double[:] c = np.zeros(N1 + N2 + 2) # coincidences
cdef double[:] mp = np.ones(N1 + N2 + 2) # multiplicity
+ cdef double interval = t_end - t_start
cdef double tau
while i + j < N1 + N2 - 2:
if (i < N1-1) and (j == N2-1 or spikes1[i+1] < spikes2[j+1]):
i += 1
n += 1
- tau = get_tau(spikes1, spikes2, i, j, max_tau)
+ tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
st[n] = spikes1[i]
if j > -1 and spikes1[i]-spikes2[j] < tau:
# coincidence between the current spike and the previous spike
@@ -392,7 +393,7 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2,
elif (j < N2-1) and (i == N1-1 or spikes1[i+1] > spikes2[j+1]):
j += 1
n += 1
- tau = get_tau(spikes1, spikes2, i, j, max_tau)
+ tau = get_tau(spikes1, spikes2, i, j, interval, max_tau)
st[n] = spikes2[j]
if i > -1 and spikes2[j]-spikes1[i] < tau:
# coincidence between the current spike and the previous spike
@@ -415,9 +416,13 @@ def coincidence_cython(double[:] spikes1, double[:] spikes2,
st[0] = t_start
st[len(st)-1] = t_end
- c[0] = c[1]
- c[len(c)-1] = c[len(c)-2]
- mp[0] = mp[1]
- mp[len(mp)-1] = mp[len(mp)-2]
+ if N1 + N2 > 0:
+ c[0] = c[1]
+ c[len(c)-1] = c[len(c)-2]
+ mp[0] = mp[1]
+ mp[len(mp)-1] = mp[len(mp)-2]
+ else:
+ c[0] = 1
+ c[1] = 1
return st, c, mp
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index 1fd8c42..69a420f 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -340,7 +340,7 @@ def cumulative_sync_python(spikes1, spikes2):
def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
def get_tau(spikes1, spikes2, i, j, max_tau):
- m = 1E100 # some huge number
+ m = t_end - t_start # use interval as initial tau
if i < len(spikes1)-1 and i > -1:
m = min(m, spikes1[i+1]-spikes1[i])
if j < len(spikes2)-1 and j > -1:
@@ -399,10 +399,14 @@ def coincidence_python(spikes1, spikes2, t_start, t_end, max_tau):
st[0] = t_start
st[len(st)-1] = t_end
- c[0] = c[1]
- c[len(c)-1] = c[len(c)-2]
- mp[0] = mp[1]
- mp[len(mp)-1] = mp[len(mp)-2]
+ if N1 + N2 > 0:
+ c[0] = c[1]
+ c[len(c)-1] = c[len(c)-2]
+ mp[0] = mp[1]
+ mp[len(mp)-1] = mp[len(mp)-2]
+ else:
+ c[0] = 1
+ c[1] = 1
return st, c, mp
diff --git a/pyspike/generic.py b/pyspike/generic.py
index 4f278d2..41affcb 100644
--- a/pyspike/generic.py
+++ b/pyspike/generic.py
@@ -31,6 +31,69 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
Returns:
- The averaged multi-variate distance of all pairs
"""
+
+ def divide_and_conquer(pairs1, pairs2):
+ """ recursive calls by splitting the two lists in half.
+ """
+ L1 = len(pairs1)
+ if L1 > 1:
+ dist_prof1 = divide_and_conquer(pairs1[:L1/2], pairs1[L1/2:])
+ else:
+ dist_prof1 = pair_distance_func(spike_trains[pairs1[0][0]],
+ spike_trains[pairs1[0][1]])
+ L2 = len(pairs2)
+ if L2 > 1:
+ dist_prof2 = divide_and_conquer(pairs2[:L2/2], pairs2[L2/2:])
+ else:
+ dist_prof2 = pair_distance_func(spike_trains[pairs2[0][0]],
+ spike_trains[pairs2[0][1]])
+ dist_prof1.add(dist_prof2)
+ return dist_prof1
+
+ if indices is None:
+ indices = np.arange(len(spike_trains))
+ indices = np.array(indices)
+ # check validity of indices
+ assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
+ "Invalid index list."
+ # generate a list of possible index pairs
+ pairs = [(indices[i], j) for i in range(len(indices))
+ for j in indices[i+1:]]
+
+ L = len(pairs)
+ if L > 1:
+ # recursive iteration through the list of pairs to get average profile
+ avrg_dist = divide_and_conquer(pairs[:len(pairs)/2],
+ pairs[len(pairs)/2:])
+ else:
+ avrg_dist = pair_distance_func(spike_trains[pairs[0][0]],
+ spike_trains[pairs[0][1]])
+
+ return avrg_dist, L
+
+
+############################################################
+# _generic_distance_multi
+############################################################
+def _generic_distance_multi(spike_trains, pair_distance_func,
+ indices=None, interval=None):
+ """ Internal implementation detail, don't call this function directly,
+ use isi_distance_multi or spike_distance_multi instead.
+
+ Computes the multi-variate distance for a set of spike-trains using the
+ pair_dist_func to compute pair-wise distances. That is it computes the
+ average distance of all pairs of spike-trains:
+ :math:`S(t) = 2/((N(N-1)) sum_{<i,j>} D_{i,j}`,
+ where the sum goes over all pairs <i,j>.
+ Args:
+ - spike_trains: list of spike trains
+ - pair_distance_func: function computing the distance of two spike trains
+ - indices: list of indices defining which spike trains to use,
+ if None all given spike trains are used (default=None)
+ Returns:
+ - The averaged multi-variate distance of all pairs
+ """
+
if indices is None:
indices = np.arange(len(spike_trains))
indices = np.array(indices)
@@ -40,13 +103,13 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
# generate a list of possible index pairs
pairs = [(indices[i], j) for i in range(len(indices))
for j in indices[i+1:]]
- # start with first pair
- (i, j) = pairs[0]
- average_dist = pair_distance_func(spike_trains[i], spike_trains[j])
- for (i, j) in pairs[1:]:
- current_dist = pair_distance_func(spike_trains[i], spike_trains[j])
- average_dist.add(current_dist) # add to the average
- return average_dist, len(pairs)
+
+ avrg_dist = 0.0
+ for (i, j) in pairs:
+ avrg_dist += pair_distance_func(spike_trains[i], spike_trains[j],
+ interval)
+
+ return avrg_dist/len(pairs)
############################################################
diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py
index aeab0df..5ea555d 100644
--- a/pyspike/isi_distance.py
+++ b/pyspike/isi_distance.py
@@ -3,7 +3,8 @@
# Distributed under the BSD License
from pyspike import PieceWiseConstFunc
-from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
+from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \
+ _generic_distance_matrix
############################################################
@@ -24,24 +25,25 @@ def isi_profile(spike_train1, spike_train2):
"""
# check whether the spike trains are defined for the same interval
assert spike_train1.t_start == spike_train2.t_start, \
- "Given spike trains seems not to have auxiliary spikes!"
+ "Given spike trains are not defined on the same interval!"
assert spike_train1.t_end == spike_train2.t_end, \
- "Given spike trains seems not to have auxiliary spikes!"
+ "Given spike trains are not defined on the same interval!"
# load cython implementation
try:
- from cython.cython_distance import isi_distance_cython \
- as isi_distance_impl
+ from cython.cython_profiles import isi_profile_cython \
+ as isi_profile_impl
except ImportError:
print("Warning: isi_distance_cython not found. Make sure that PySpike \
is installed by running\n 'python setup.py build_ext --inplace'!\n \
Falling back to slow python backend.")
# use python backend
from cython.python_backend import isi_distance_python \
- as isi_distance_impl
+ as isi_profile_impl
- times, values = isi_distance_impl(spike_train1.spikes, spike_train2.spikes,
- spike_train1.t_start, spike_train1.t_end)
+ times, values = isi_profile_impl(spike_train1.get_spikes_non_empty(),
+ spike_train2.get_spikes_non_empty(),
+ spike_train1.t_start, spike_train1.t_end)
return PieceWiseConstFunc(times, values)
@@ -65,7 +67,23 @@ def isi_distance(spike_train1, spike_train2, interval=None):
:returns: The isi-distance :math:`D_I`.
:rtype: double
"""
- return isi_profile(spike_train1, spike_train2).avrg(interval)
+
+ if interval is None:
+ # distance over the whole interval is requested: use specific function
+ # for optimal performance
+ try:
+ from cython.cython_distances import isi_distance_cython \
+ as isi_distance_impl
+
+ return isi_distance_impl(spike_train1.get_spikes_non_empty(),
+ spike_train2.get_spikes_non_empty(),
+ spike_train1.t_start, spike_train1.t_end)
+ except ImportError:
+ # Cython backend not available: fall back to profile averaging
+ return isi_profile(spike_train1, spike_train2).avrg(interval)
+ else:
+ # some specific interval is provided: use profile
+ return isi_profile(spike_train1, spike_train2).avrg(interval)
############################################################
@@ -112,7 +130,8 @@ def isi_distance_multi(spike_trains, indices=None, interval=None):
:returns: The time-averaged multivariate ISI distance :math:`D_I`
:rtype: double
"""
- return isi_profile_multi(spike_trains, indices).avrg(interval)
+ return _generic_distance_multi(spike_trains, isi_distance, indices,
+ interval)
############################################################
diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py
index cc620d4..ac2d260 100644
--- a/pyspike/spike_distance.py
+++ b/pyspike/spike_distance.py
@@ -3,7 +3,8 @@
# Distributed under the BSD License
from pyspike import PieceWiseLinFunc
-from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
+from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \
+ _generic_distance_matrix
############################################################
@@ -24,26 +25,27 @@ def spike_profile(spike_train1, spike_train2):
"""
# check whether the spike trains are defined for the same interval
assert spike_train1.t_start == spike_train2.t_start, \
- "Given spike trains seems not to have auxiliary spikes!"
+ "Given spike trains are not defined on the same interval!"
assert spike_train1.t_end == spike_train2.t_end, \
- "Given spike trains seems not to have auxiliary spikes!"
+ "Given spike trains are not defined on the same interval!"
# cython implementation
try:
- from cython.cython_distance import spike_distance_cython \
- as spike_distance_impl
+ from cython.cython_profiles import spike_profile_cython \
+ as spike_profile_impl
except ImportError:
- print("Warning: spike_distance_cython not found. Make sure that \
+ print("Warning: spike_profile_cython not found. Make sure that \
PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \
Falling back to slow python backend.")
# use python backend
from cython.python_backend import spike_distance_python \
- as spike_distance_impl
+ as spike_profile_impl
+
+ times, y_starts, y_ends = spike_profile_impl(
+ spike_train1.get_spikes_non_empty(),
+ spike_train2.get_spikes_non_empty(),
+ spike_train1.t_start, spike_train1.t_end)
- times, y_starts, y_ends = spike_distance_impl(spike_train1.spikes,
- spike_train2.spikes,
- spike_train1.t_start,
- spike_train1.t_end)
return PieceWiseLinFunc(times, y_starts, y_ends)
@@ -67,7 +69,22 @@ def spike_distance(spike_train1, spike_train2, interval=None):
:rtype: double
"""
- return spike_profile(spike_train1, spike_train2).avrg(interval)
+ if interval is None:
+ # distance over the whole interval is requested: use specific function
+ # for optimal performance
+ try:
+ from cython.cython_distances import spike_distance_cython \
+ as spike_distance_impl
+ return spike_distance_impl(spike_train1.get_spikes_non_empty(),
+ spike_train2.get_spikes_non_empty(),
+ spike_train1.t_start,
+ spike_train1.t_end)
+ except ImportError:
+ # Cython backend not available: fall back to average profile
+ return spike_profile(spike_train1, spike_train2).avrg(interval)
+ else:
+ # some specific interval is provided: compute the whole profile
+ return spike_profile(spike_train1, spike_train2).avrg(interval)
############################################################
@@ -117,7 +134,8 @@ def spike_distance_multi(spike_trains, indices=None, interval=None):
:returns: The averaged multi-variate spike distance :math:`D_S`.
:rtype: double
"""
- return spike_profile_multi(spike_trains, indices).avrg(interval)
+ return _generic_distance_multi(spike_trains, spike_distance, indices,
+ interval)
############################################################
diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py
index 9d2e363..40d98d2 100644
--- a/pyspike/spike_sync.py
+++ b/pyspike/spike_sync.py
@@ -3,6 +3,7 @@
# Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
# Distributed under the BSD License
+import numpy as np
from functools import partial
from pyspike import DiscreteFunc
from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
@@ -29,35 +30,68 @@ def spike_sync_profile(spike_train1, spike_train2, max_tau=None):
"""
# check whether the spike trains are defined for the same interval
assert spike_train1.t_start == spike_train2.t_start, \
- "Given spike trains seems not to have auxiliary spikes!"
+ "Given spike trains are not defined on the same interval!"
assert spike_train1.t_end == spike_train2.t_end, \
- "Given spike trains seems not to have auxiliary spikes!"
+ "Given spike trains are not defined on the same interval!"
# cython implementation
try:
- from cython.cython_distance import coincidence_cython \
- as coincidence_impl
+ from cython.cython_profiles import coincidence_profile_cython \
+ as coincidence_profile_impl
except ImportError:
print("Warning: spike_distance_cython not found. Make sure that \
PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \
Falling back to slow python backend.")
# use python backend
from cython.python_backend import coincidence_python \
- as coincidence_impl
+ as coincidence_profile_impl
if max_tau is None:
max_tau = 0.0
- times, coincidences, multiplicity = coincidence_impl(spike_train1.spikes,
- spike_train2.spikes,
- spike_train1.t_start,
- spike_train1.t_end,
- max_tau)
+ times, coincidences, multiplicity \
+ = coincidence_profile_impl(spike_train1.spikes, spike_train2.spikes,
+ spike_train1.t_start, spike_train1.t_end,
+ max_tau)
return DiscreteFunc(times, coincidences, multiplicity)
############################################################
+# _spike_sync_values
+############################################################
+def _spike_sync_values(spike_train1, spike_train2, interval, max_tau):
+ """" Internal function. Computes the summed coincidences and multiplicity
+ for spike synchronization of the two given spike trains.
+
+ Do not call this function directly, use `spike_sync` or `spike_sync_multi`
+ instead.
+ """
+ if interval is None:
+ # distance over the whole interval is requested: use specific function
+ # for optimal performance
+ try:
+ from cython.cython_distances import coincidence_value_cython \
+ as coincidence_value_impl
+ if max_tau is None:
+ max_tau = 0.0
+ c, mp = coincidence_value_impl(spike_train1.spikes,
+ spike_train2.spikes,
+ spike_train1.t_start,
+ spike_train1.t_end,
+ max_tau)
+ return c, mp
+ except ImportError:
+ # Cython backend not available: fall back to profile averaging
+ return spike_sync_profile(spike_train1, spike_train2,
+ max_tau).integral(interval)
+ else:
+ # some specific interval is provided: use profile
+ return spike_sync_profile(spike_train1, spike_train2,
+ max_tau).integral(interval)
+
+
+############################################################
# spike_sync
############################################################
def spike_sync(spike_train1, spike_train2, interval=None, max_tau=None):
@@ -80,8 +114,8 @@ def spike_sync(spike_train1, spike_train2, interval=None, max_tau=None):
:rtype: `double`
"""
- return spike_sync_profile(spike_train1, spike_train2,
- max_tau).avrg(interval)
+ c, mp = _spike_sync_values(spike_train1, spike_train2, interval, max_tau)
+ return 1.0*c/mp
############################################################
@@ -131,8 +165,25 @@ def spike_sync_multi(spike_trains, indices=None, interval=None, max_tau=None):
:rtype: double
"""
- return spike_sync_profile_multi(spike_trains, indices,
- max_tau).avrg(interval)
+ if indices is None:
+ indices = np.arange(len(spike_trains))
+ indices = np.array(indices)
+ # check validity of indices
+ assert (indices < len(spike_trains)).all() and (indices >= 0).all(), \
+ "Invalid index list."
+ # generate a list of possible index pairs
+ pairs = [(indices[i], j) for i in range(len(indices))
+ for j in indices[i+1:]]
+
+ coincidence = 0.0
+ mp = 0.0
+ for (i, j) in pairs:
+ c, m = _spike_sync_values(spike_trains[i], spike_trains[j],
+ interval, max_tau)
+ coincidence += c
+ mp += m
+
+ return coincidence/mp
############################################################
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-29 18:02:33 +0000