summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMario Mulansky <mario.mulansky@gmx.net>2015-12-22 18:15:59 +0100
committerMario Mulansky <mario.mulansky@gmx.net>2015-12-22 18:15:59 +0100
commite32272f4540de347abcc548a94239b625458b3a6 (patch)
tree409b60357e7dae2ff40d0e96cb9345d5673d431f
parent94c5fd007d33a38f3c9d1121749cb6ffb162394c (diff)
changed edge correction for single spikes
Spike trains with single spikes now only get auxiliary spikes at the edges for the SPIKE distance instead of real spikes before.
-rw-r--r--pyspike/SpikeTrain.py2
-rw-r--r--pyspike/cython/cython_distances.pyx89
-rw-r--r--pyspike/cython/cython_profiles.pyx67
-rw-r--r--pyspike/cython/python_backend.py68
-rw-r--r--test/test_empty.py16
5 files changed, 130 insertions, 112 deletions
diff --git a/pyspike/SpikeTrain.py b/pyspike/SpikeTrain.py
index 4b59a5d..19f2419 100644
--- a/pyspike/SpikeTrain.py
+++ b/pyspike/SpikeTrain.py
@@ -68,7 +68,7 @@ class SpikeTrain(object):
"""Returns the spikes of this spike train with auxiliary spikes in case
of empty spike trains.
"""
- if len(self.spikes) < 2:
+ if len(self.spikes) < 1:
return np.unique(np.insert([self.t_start, self.t_end], 1,
self.spikes))
else:
diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx
index 6c6e7e5..7dc9cb9 100644
--- a/pyspike/cython/cython_distances.pyx
+++ b/pyspike/cython/cython_distances.pyx
@@ -55,20 +55,27 @@ def isi_distance_cython(double[:] s1, double[:] s2,
N2 = len(s2)
# first interspike interval - check if a spike exists at the start time
+ # and also account for spike trains with single spikes
if s1[0] > t_start:
- # edge correction
- nu1 = fmax(s1[0]-t_start, s1[1]-s1[0])
+ # edge correction for the first interspike interval:
+ # take the maximum of the distance from the beginning to the first
+ # spike and the interval between the first two spikes.
+ # if there is only one spike, take the its distance to the beginning
+ nu1 = fmax(s1[0]-t_start, s1[1]-s1[0]) if N1 > 1 else s1[0]-t_start
index1 = -1
else:
- nu1 = s1[1]-s1[0]
+ # if the first spike is exactly at the start, take the distance
+ # to the next spike. If this is the only spike, take the distance to
+ # the end.
+ nu1 = s1[1]-s1[0] if N1 > 1 else t_end-s1[0]
index1 = 0
if s2[0] > t_start:
- # edge correction
- nu2 = fmax(s2[0]-t_start, s2[1]-s2[0])
+ # edge correction as above
+ nu2 = fmax(s2[0]-t_start, s2[1]-s2[0]) if N2 > 1 else s2[0]-t_start
index2 = -1
else:
- nu2 = s2[1]-s2[0]
+ nu2 = s2[1]-s2[0] if N2 > 1 else t_end-s2[0]
index2 = 0
last_t = t_start
@@ -86,8 +93,12 @@ def isi_distance_cython(double[:] s1, double[:] s2,
if index1 < N1-1:
nu1 = s1[index1+1]-s1[index1]
else:
- # edge correction
- nu1 = fmax(t_end-s1[index1], nu1)
+ # edge correction for the last ISI:
+ # take the max of the distance of the last
+ # spike to the end and the previous ISI. If there was only
+ # one spike, always take the distance to the end.
+ nu1 = fmax(t_end-s1[index1], nu1) if N1 > 1 \
+ else t_end-s1[index1]
elif (index2 < N2-1) and ((index1 == N1-1) or
(s1[index1+1] > s2[index2+1])):
index2 += 1
@@ -95,8 +106,9 @@ def isi_distance_cython(double[:] s1, double[:] s2,
if index2 < N2-1:
nu2 = s2[index2+1]-s2[index2]
else:
- # edge correction
- nu2 = fmax(t_end-s2[index2], nu2)
+ # edge correction for the end as above
+ nu2 = fmax(t_end-s2[index2], nu2) if N2 > 1 \
+ else t_end-s2[index2]
else: # s1[index1+1] == s2[index2+1]
index1 += 1
index2 += 1
@@ -104,13 +116,15 @@ def isi_distance_cython(double[:] s1, double[:] s2,
if index1 < N1-1:
nu1 = s1[index1+1]-s1[index1]
else:
- # edge correction
- nu1 = fmax(t_end-s1[index1], nu1)
+ # edge correction for the end as above
+ nu1 = fmax(t_end-s1[index1], nu1) if N1 > 1 \
+ else t_end-s1[index1]
if index2 < N2-1:
nu2 = s2[index2+1]-s2[index2]
else:
- # edge correction
- nu2 = fmax(t_end-s2[index2], nu2)
+ # edge correction for the end as above
+ nu2 = fmax(t_end-s2[index2], nu2) if N2 > 1 \
+ else t_end-s2[index2]
# compute the corresponding isi-distance
isi_value += curr_isi * (curr_t - last_t)
curr_isi = fabs(nu1 - nu2) / fmax(nu1, nu2)
@@ -184,16 +198,18 @@ def spike_distance_cython(double[:] t1, double[:] t2,
N1 = len(t1)
N2 = len(t2)
- # with nogil: # release the interpreter to allow multithreading
- if True:
+ # we can assume at least one spikes per spike train
+ assert N1 > 0
+ assert N2 > 0
+
+
+ with nogil: # release the interpreter to allow multithreading
t_last = t_start
- # t_p1 = t_start
- # t_p2 = t_start
# auxiliary spikes for edge correction - consistent with first/last ISI
- t_aux1[0] = fmin(t_start, t1[0]-(t1[1]-t1[0]))
- t_aux1[1] = fmax(t_end, t1[N1-1]+(t1[N1-1]-t1[N1-2]))
- t_aux2[0] = fmin(t_start, t2[0]-(t2[1]-t2[0]))
- t_aux2[1] = fmax(t_end, t2[N2-1]+(t2[N2-1]-t2[N2-2]))
+ t_aux1[0] = fmin(t_start, 2*t1[0]-t1[1]) if N1 > 1 else t_start
+ t_aux1[1] = fmax(t_end, 2*t1[N1-1]-t1[N1-2]) if N1 > 1 else t_end
+ t_aux2[0] = fmin(t_start, 2*t2[0]-t2[1]) if N2 > 1 else t_start
+ t_aux2[1] = fmax(t_end, 2*t2[N2-1]+-t2[N2-2]) if N2 > 1 else t_end
# print "aux spikes %.15f, %.15f ; %.15f, %.15f" % (t_aux1[0], t_aux1[1], t_aux2[0], t_aux2[1])
t_p1 = t_start if (t1[0] == t_start) else t_aux1[0]
t_p2 = t_start if (t2[0] == t_start) else t_aux2[0]
@@ -201,17 +217,16 @@ def spike_distance_cython(double[:] t1, double[:] t2,
# dt_p1 = t2[0]-t_start
t_f1 = t1[0]
dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_aux2[0], t_aux2[1])
- isi1 = fmax(t_f1-t_start, t1[1]-t1[0])
+ isi1 = fmax(t_f1-t_start, t1[1]-t1[0]) if N1 > 1 else t_f1-t_start
dt_p1 = dt_f1
# s1 = dt_p1*(t_f1-t_start)/isi1
s1 = dt_p1
index1 = -1
- else:
- t_f1 = t1[1]
+ else: # t1[0] == t_start
+ t_f1 = t1[1] if N1 > 1 else t_end
dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_aux2[0], t_aux2[1])
- # dt_p1 = t_start-t_p2 # 0.0
dt_p1 = get_min_dist_cython(t_p1, t2, N2, 0, t_aux2[0], t_aux2[1])
- isi1 = t1[1]-t1[0]
+ isi1 = t_f1-t1[0]
s1 = dt_p1
index1 = 0
if t2[0] > t_start:
@@ -219,16 +234,16 @@ def spike_distance_cython(double[:] t1, double[:] t2,
t_f2 = t2[0]
dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1])
dt_p2 = dt_f2
- isi2 = fmax(t_f2-t_start, t2[1]-t2[0])
+ isi2 = fmax(t_f2-t_start, t2[1]-t2[0]) if N2 > 1 else t_f2-t_start
# s2 = dt_p2*(t_f2-t_start)/isi2
s2 = dt_p2
index2 = -1
- else:
- t_f2 = t2[1]
+ else: # t2[0] == t_start
+ t_f2 = t2[1] if N2 > 1 else t_end
dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1])
# dt_p2 = t_start-t_p1 # 0.0
dt_p2 = get_min_dist_cython(t_p2, t1, N1, 0, t_aux1[0], t_aux1[1])
- isi2 = t2[1]-t2[0]
+ isi2 = t_f2-t2[0]
s2 = dt_p2
index2 = 0
@@ -263,7 +278,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
s1 = dt_p1
else:
dt_f1 = dt_p1
- isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) if N1 > 1 \
+ else t_end-t1[N1-1]
# s1 needs adjustment due to change of isi1
# s1 = dt_p1*(t_end-t1[N1-1])/isi1
# Eero's correction: no adjustment
@@ -296,7 +312,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
s2 = dt_p2
else:
dt_f2 = dt_p2
- isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) if N2 > 1 \
+ else t_end-t2[N2-1]
# s2 needs adjustment due to change of isi2
# s2 = dt_p2*(t_end-t2[N2-1])/isi2
# Eero's correction: no adjustment
@@ -322,7 +339,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
else:
t_f1 = t_aux1[1]
dt_f1 = dt_p1
- isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) if N1 > 1 \
+ else t_end-t1[N1-1]
if index2 < N2-1:
t_f2 = t2[index2+1]
dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
@@ -331,7 +349,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
else:
t_f2 = t_aux2[1]
dt_f2 = dt_p2
- isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) if N2 > 1 \
+ else t_end-t2[N2-1]
index += 1
t_last = t_curr
# isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx
index f08de6e..4a42cdb 100644
--- a/pyspike/cython/cython_profiles.pyx
+++ b/pyspike/cython/cython_profiles.pyx
@@ -63,18 +63,18 @@ def isi_profile_cython(double[:] s1, double[:] 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])
+ nu1 = fmax(s1[0]-t_start, s1[1]-s1[0]) if N1 > 1 else s1[0]-t_start
index1 = -1
else:
- nu1 = s1[1]-s1[0]
+ nu1 = s1[1]-s1[0] if N1 > 1 else t_end-s1[0]
index1 = 0
if s2[0] > t_start:
# edge correction
- nu2 = fmax(s2[0]-t_start, s2[1]-s2[0])
+ nu2 = fmax(s2[0]-t_start, s2[1]-s2[0]) if N2 > 1 else s2[0]-t_start
index2 = -1
else:
- nu2 = s2[1]-s2[0]
+ nu2 = s2[1]-s2[0] if N2 > 1 else t_end-s2[0]
index2 = 0
isi_values[0] = fabs(nu1-nu2)/fmax(nu1, nu2)
@@ -92,7 +92,8 @@ def isi_profile_cython(double[:] s1, double[:] s2,
nu1 = s1[index1+1]-s1[index1]
else:
# edge correction
- nu1 = fmax(t_end-s1[index1], nu1)
+ nu1 = fmax(t_end-s1[index1], nu1) if N1 > 1 \
+ else t_end-s1[index1]
elif (index2 < N2-1) and ((index1 == N1-1) or
(s1[index1+1] > s2[index2+1])):
index2 += 1
@@ -101,7 +102,8 @@ def isi_profile_cython(double[:] s1, double[:] s2,
nu2 = s2[index2+1]-s2[index2]
else:
# edge correction
- nu2 = fmax(t_end-s2[index2], nu2)
+ nu2 = fmax(t_end-s2[index2], nu2) if N2 > 1 \
+ else t_end-s2[index2]
else: # s1[index1+1] == s2[index2+1]
index1 += 1
index2 += 1
@@ -110,12 +112,14 @@ def isi_profile_cython(double[:] s1, double[:] s2,
nu1 = s1[index1+1]-s1[index1]
else:
# edge correction
- nu1 = fmax(t_end-s1[index1], nu1)
+ nu1 = fmax(t_end-s1[index1], nu1) if N1 > 1 \
+ else t_end-s1[index1]
if index2 < N2-1:
nu2 = s2[index2+1]-s2[index2]
else:
# edge correction
- nu2 = fmax(t_end-s2[index2], nu2)
+ nu2 = fmax(t_end-s2[index2], nu2) if N2 > 1 \
+ else t_end-s2[index2]
# compute the corresponding isi-distance
isi_values[index] = fabs(nu1 - nu2) / fmax(nu1, nu2)
index += 1
@@ -191,9 +195,9 @@ def spike_profile_cython(double[:] t1, double[:] t2,
N1 = len(t1)
N2 = len(t2)
- # we can assume at least two spikes per spike train
- assert N1 > 1
- assert N2 > 1
+ # we can assume at least one spikes per spike train
+ assert N1 > 0
+ assert N2 > 0
spike_events = np.empty(N1+N2+2)
@@ -205,26 +209,26 @@ def spike_profile_cython(double[:] t1, double[:] t2,
# t_p1 = t_start
# t_p2 = t_start
# auxiliary spikes for edge correction - consistent with first/last ISI
- t_aux1[0] = fmin(t_start, t1[0]-(t1[1]-t1[0]))
- t_aux1[1] = fmax(t_end, t1[N1-1]+(t1[N1-1]-t1[N1-2]))
- t_aux2[0] = fmin(t_start, t2[0]-(t2[1]-t2[0]))
- t_aux2[1] = fmax(t_end, t2[N2-1]+(t2[N2-1]-t2[N2-2]))
+ t_aux1[0] = fmin(t_start, 2*t1[0]-t1[1]) if N1 > 1 else t_start
+ t_aux1[1] = fmax(t_end, 2*t1[N1-1]-t1[N1-2]) if N1 > 1 else t_end
+ t_aux2[0] = fmin(t_start, 2*t2[0]-t2[1]) if N2 > 1 else t_start
+ t_aux2[1] = fmax(t_end, 2*t2[N2-1]-t2[N2-2]) if N2 > 1 else t_end
t_p1 = t_start if (t1[0] == t_start) else t_aux1[0]
t_p2 = t_start if (t2[0] == t_start) else t_aux2[0]
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_aux2[0], t_aux2[1])
- isi1 = fmax(t_f1-t_start, t1[1]-t1[0])
+ isi1 = fmax(t_f1-t_start, t1[1]-t1[0]) if N1 > 1 else t_f1-t_start
dt_p1 = dt_f1
# s1 = dt_p1*(t_f1-t_start)/isi1
s1 = dt_p1
index1 = -1
else:
- t_f1 = t1[1]
+ t_f1 = t1[1] if N1 > 1 else t_end
dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_aux2[0], t_aux2[1])
- dt_p1 = 0.0
- isi1 = t1[1]-t1[0]
+ dt_p1 = get_min_dist_cython(t_p1, t2, N2, 0, t_aux2[0], t_aux2[1])
+ isi1 = t_f1-t1[0]
s1 = dt_p1
index1 = 0
if t2[0] > t_start:
@@ -232,15 +236,15 @@ def spike_profile_cython(double[:] t1, double[:] t2,
t_f2 = t2[0]
dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1])
dt_p2 = dt_f2
- isi2 = fmax(t_f2-t_start, t2[1]-t2[0])
+ isi2 = fmax(t_f2-t_start, t2[1]-t2[0]) if N2 > 1 else t_f2-t_start
# s2 = dt_p2*(t_f2-t_start)/isi2
s2 = dt_p2
index2 = -1
else:
- t_f2 = t2[1]
+ t_f2 = t2[1] if N2 > 1 else t_end
dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_aux1[0], t_aux1[1])
- dt_p2 = 0.0
- isi2 = t2[1]-t2[0]
+ dt_p2 = get_min_dist_cython(t_p2, t1, N1, 0, t_aux1[0], t_aux1[1])
+ isi2 = t_f2-t2[0]
s2 = dt_p2
index2 = 0
@@ -273,7 +277,8 @@ def spike_profile_cython(double[:] t1, double[:] t2,
s1 = dt_p1
else:
dt_f1 = dt_p1
- isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) if N1 > 1 \
+ else t_end-t1[N1-1]
# s1 needs adjustment due to change of isi1
# s1 = dt_p1*(t_end-t1[N1-1])/isi1
# Eero's correction: no adjustment
@@ -305,7 +310,8 @@ def spike_profile_cython(double[:] t1, double[:] t2,
s2 = dt_p2
else:
dt_f2 = dt_p2
- isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) if N2 > 1 \
+ else t_end-t2[N2-1]
# s2 needs adjustment due to change of isi2
# s2 = dt_p2*(t_end-t2[N2-1])/isi2
# Eero's correction: no adjustment
@@ -330,7 +336,8 @@ def spike_profile_cython(double[:] t1, double[:] t2,
else:
t_f1 = t_aux1[1]
dt_f1 = dt_p1
- isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ isi1 = fmax(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) if N1 > 1 \
+ else t_end-t1[N1-1]
if index2 < N2-1:
t_f2 = t2[index2+1]
dt_f2 = get_min_dist_cython(t_f2, t1, N1, index1,
@@ -339,18 +346,14 @@ def spike_profile_cython(double[:] t1, double[:] t2,
else:
t_f2 = t_aux2[1]
dt_f2 = dt_p2
- isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) if N2 > 1 \
+ else t_end-t2[N2-1]
index += 1
# the last event is the interval end
if spike_events[index-1] == t_end:
index -= 1
else:
spike_events[index] = t_end
- # the ending value of the last interval
- 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
s1 = dt_f1
s2 = dt_f2
y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1, isi2)
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index cf898d7..a1e3a65 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -28,17 +28,17 @@ def isi_distance_python(s1, s2, t_start, t_end):
isi_values = np.empty(len(spike_events) - 1)
if s1[0] > t_start:
# edge correction
- nu1 = max(s1[0] - t_start, s1[1] - s1[0])
+ nu1 = max(s1[0] - t_start, s1[1] - s1[0]) if N1 > 1 else s1[0]-t_start
index1 = -1
else:
- nu1 = s1[1] - s1[0]
+ nu1 = s1[1] - s1[0] if N1 > 1 else t_end-s1[0]
index1 = 0
if s2[0] > t_start:
# edge correction
- nu2 = max(s2[0] - t_start, s2[1] - s2[0])
+ nu2 = max(s2[0] - t_start, s2[1] - s2[0]) if N2 > 1 else s2[0]-t_start
index2 = -1
else:
- nu2 = s2[1] - s2[0]
+ nu2 = s2[1] - s2[0] if N2 > 1 else t_end-s2[0]
index2 = 0
isi_values[0] = abs(nu1 - nu2) / max(nu1, nu2)
@@ -52,7 +52,8 @@ def isi_distance_python(s1, s2, t_start, t_end):
nu1 = s1[index1+1]-s1[index1]
else:
# edge correction
- nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2])
+ nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2]) if N1 > 1 \
+ else t_end-s1[N1-1]
elif (index2 < N2-1) and (index1 == N1-1 or
s1[index1+1] > s2[index2+1]):
@@ -62,7 +63,8 @@ def isi_distance_python(s1, s2, t_start, t_end):
nu2 = s2[index2+1]-s2[index2]
else:
# edge correction
- nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2])
+ nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2]) if N2 > 1 \
+ else t_end-s2[N2-1]
else: # s1[index1 + 1] == s2[index2 + 1]
index1 += 1
@@ -72,12 +74,14 @@ def isi_distance_python(s1, s2, t_start, t_end):
nu1 = s1[index1+1]-s1[index1]
else:
# edge correction
- nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2])
+ nu1 = max(t_end-s1[N1-1], s1[N1-1]-s1[N1-2]) if N1 > 1 \
+ else t_end-s1[N1-1]
if index2 < N2-1:
nu2 = s2[index2+1]-s2[index2]
else:
# edge correction
- nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2])
+ nu2 = max(t_end-s2[N2-1], s2[N2-1]-s2[N2-2]) if N2 > 1 \
+ else t_end-s2[N2-1]
# compute the corresponding isi-distance
isi_values[index] = abs(nu1 - nu2) / \
max(nu1, nu2)
@@ -146,10 +150,10 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
t_aux1 = np.zeros(2)
t_aux2 = np.zeros(2)
- t_aux1[0] = min(t_start, t1[0]-(t1[1]-t1[0]))
- t_aux1[1] = max(t_end, t1[N1-1]+(t1[N1-1]-t1[N1-2]))
- t_aux2[0] = min(t_start, t2[0]-(t2[1]-t2[0]))
- t_aux2[1] = max(t_end, t2[N2-1]+(t2[N2-1]-t2[N2-2]))
+ t_aux1[0] = min(t_start, t1[0]-(t1[1]-t1[0])) if N1 > 1 else t_start
+ t_aux1[1] = max(t_end, t1[N1-1]+(t1[N1-1]-t1[N1-2])) if N1 > 1 else t_end
+ t_aux2[0] = min(t_start, t2[0]-(t2[1]-t2[0])) if N2 > 1 else t_start
+ t_aux2[1] = max(t_end, t2[N2-1]+(t2[N2-1]-t2[N2-2])) if N2 > 1 else t_end
t_p1 = t_start if (t1[0] == t_start) else t_aux1[0]
t_p2 = t_start if (t2[0] == t_start) else t_aux2[0]
@@ -160,16 +164,16 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
t_f1 = t1[0]
dt_f1 = get_min_dist(t_f1, t2, 0, t_aux2[0], t_aux2[1])
dt_p1 = dt_f1
- isi1 = max(t_f1-t_start, t1[1]-t1[0])
+ isi1 = max(t_f1-t_start, t1[1]-t1[0]) if N1 > 1 else t_f1-t_start
# s1 = dt_p1*(t_f1-t_start)/isi1
s1 = dt_p1
index1 = -1
else:
# dt_p1 = t_start-t_p2
+ t_f1 = t1[1] if N1 > 1 else t_end
dt_p1 = get_min_dist(t_p1, t2, 0, t_aux2[0], t_aux2[1])
- t_f1 = t1[1]
dt_f1 = get_min_dist(t_f1, t2, 0, t_aux2[0], t_aux2[1])
- isi1 = t1[1]-t1[0]
+ isi1 = t_f1-t1[0]
s1 = dt_p1
index1 = 0
if t2[0] > t_start:
@@ -177,23 +181,18 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
t_f2 = t2[0]
dt_f2 = get_min_dist(t_f2, t1, 0, t_aux1[0], t_aux1[1])
dt_p2 = dt_f2
- isi2 = max(t_f2-t_start, t2[1]-t2[0])
+ isi2 = max(t_f2-t_start, t2[1]-t2[0]) if N2 > 1 else t_f2-t_start
# s2 = dt_p2*(t_f2-t_start)/isi2
s2 = dt_p2
index2 = -1
else:
- dt_p2 = t_start-t_p1
+ t_f2 = t2[1] if N2 > 1 else t_end
dt_p2 = get_min_dist(t_p2, t1, 0, t_aux1[0], t_aux1[1])
- t_f2 = t2[1]
dt_f2 = get_min_dist(t_f2, t1, 0, t_aux1[0], t_aux1[1])
- isi2 = t2[1]-t2[0]
+ isi2 = t_f2-t2[0]
s2 = dt_p2
index2 = 0
- # print "t_p1:", repr(t_p1), ", t_f1:", repr(t_f1), ", dt_p1:", repr(dt_p1), ", dt_f1:", repr(dt_f1)
- # print "t_p2:", repr(t_p2), ", t_f2:", repr(t_f2), ", dt_p2:", repr(dt_p2), ", dt_f2:", repr(dt_f2)
- # print "s1: ", repr(s1), ", s2:", repr(s2)
-
y_starts[0] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
index = 1
@@ -221,7 +220,8 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
s1 = dt_p1
else:
dt_f1 = dt_p1
- isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) if N1 > 1 \
+ else t_end-t1[N1-1]
# s1 needs adjustment due to change of isi1
# s1 = dt_p1*(t_end-t1[N1-1])/isi1
# Eero's correction: no adjustment
@@ -250,7 +250,8 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
s2 = dt_p2
else:
dt_f2 = dt_p2
- isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) if N2 > 1 \
+ else t_end-t2[N2-1]
# s2 needs adjustment due to change of isi2
# s2 = dt_p2*(t_end-t2[N2-1])/isi2
# Eero's adjustment: no correction
@@ -274,7 +275,8 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
else:
t_f1 = t_aux1[1]
dt_f1 = dt_p1
- isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2])
+ isi1 = max(t_end-t1[N1-1], t1[N1-1]-t1[N1-2]) if N1 > 1 \
+ else t_end-t1[N1-1]
if index2 < N2-1:
t_f2 = t2[index2+1]
dt_f2 = get_min_dist(t_f2, t1, index1, t_aux1[0], t_aux1[1])
@@ -282,29 +284,19 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
else:
t_f2 = t_aux2[1]
dt_f2 = dt_p2
- isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
+ isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2]) if N2 > 1 \
+ else t_end-t2[N2-1]
index += 1
- # print "t_p1:", repr(t_p1), ", t_f1:", repr(t_f1), ", dt_p1:", repr(dt_p1), ", dt_f1:", repr(dt_f1)
- # print "t_p2:", repr(t_p2), ", t_f2:", repr(t_f2), ", dt_p2:", repr(dt_p2), ", dt_f2:", repr(dt_f2)
- # print "s1: ", repr(s1), ", s2:", repr(s2)
-
# the last event is the interval end
if spike_events[index-1] == t_end:
index -= 1
else:
spike_events[index] = t_end
- # the ending value of the last interval
- 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_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
- # print "t_p1:", repr(t_p1), ", t_f1:", repr(t_f1), ", dt_p1:", repr(dt_p1), ", dt_f1:", repr(dt_f1)
- # print "t_p2:", repr(t_p2), ", t_f2:", repr(t_f2), ", dt_p2:", repr(dt_p2), ", dt_f2:", repr(dt_f2)
- # print "s1: ", repr(s1), ", s2:", repr(s2)
-
# use only the data added above
# could be less than original length due to equal spike times
return spike_events[:index+1], y_starts[:index], y_ends[:index]
diff --git a/test/test_empty.py b/test/test_empty.py
index 5a0042f..4d0a5cf 100644
--- a/test/test_empty.py
+++ b/test/test_empty.py
@@ -24,7 +24,9 @@ def test_get_non_empty():
st = SpikeTrain([0.5, ], edges=(0.0, 1.0))
spikes = st.get_spikes_non_empty()
- assert_array_equal(spikes, [0.0, 0.5, 1.0])
+ # assert_array_equal(spikes, [0.0, 0.5, 1.0])
+ # spike trains with one spike don't get edge spikes anymore
+ assert_array_equal(spikes, [0.5, ])
def test_isi_empty():
@@ -70,21 +72,23 @@ def test_spike_empty():
st1 = SpikeTrain([], edges=(0.0, 1.0))
st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0))
d = spk.spike_distance(st1, st2)
- d_expect = 0.4*0.4*1.0/(0.4+1.0)**2 + 0.6*0.4*1.0/(0.6+1.0)**2
+ d_expect = 2*0.4*0.4*1.0/(0.4+1.0)**2 + 2*0.6*0.4*1.0/(0.6+1.0)**2
assert_almost_equal(d, d_expect, decimal=15)
prof = spk.spike_profile(st1, st2)
assert_equal(d, prof.avrg())
assert_array_equal(prof.x, [0.0, 0.4, 1.0])
- assert_array_almost_equal(prof.y1, [0.0, 2*0.4*1.0/(0.6+1.0)**2],
+ assert_array_almost_equal(prof.y1, [2*0.4*1.0/(0.4+1.0)**2,
+ 2*0.4*1.0/(0.6+1.0)**2],
decimal=15)
- assert_array_almost_equal(prof.y2, [2*0.4*1.0/(0.4+1.0)**2, 0.0],
+ assert_array_almost_equal(prof.y2, [2*0.4*1.0/(0.4+1.0)**2,
+ 2*0.4*1.0/(0.6+1.0)**2],
decimal=15)
st1 = SpikeTrain([0.6, ], edges=(0.0, 1.0))
st2 = SpikeTrain([0.4, ], edges=(0.0, 1.0))
d = spk.spike_distance(st1, st2)
- s1 = np.array([0.0, 0.4*0.2/0.6, 0.2, 0.0])
- s2 = np.array([0.0, 0.2, 0.2*0.4/0.6, 0.0])
+ s1 = np.array([0.2, 0.2, 0.2, 0.2])
+ s2 = np.array([0.2, 0.2, 0.2, 0.2])
isi1 = np.array([0.6, 0.6, 0.4])
isi2 = np.array([0.4, 0.6, 0.6])
expected_y1 = (s1[:-1]*isi2+s2[:-1]*isi1) / (0.5*(isi1+isi2)**2)