summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-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)
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-28 01:10:23 +0000