summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorMario Mulansky <mario.mulansky@gmx.net>2015-12-14 17:25:53 +0100
committerMario Mulansky <mario.mulansky@gmx.net>2015-12-14 17:25:53 +0100
commitf7b90618f01d4dbf015b3d21c6c06dec8d26bd9f (patch)
tree799a2fa0e0f7558ece400de15ca8d1358565c066
parentd985f3a8de6ae840c8a127653b3d9affb1a8aa40 (diff)
parent9061f2a0c13134e53f937d730295a421fd671ea3 (diff)
Merge pull request #20 from mariomulansky/develop
Develop merge for version 0.4
-rw-r--r--.travis.yml3
-rw-r--r--Changelog6
-rw-r--r--Contributors.txt1
-rw-r--r--MANIFEST.in7
-rw-r--r--doc/conf.py4
-rw-r--r--examples/multivariate.py4
-rw-r--r--examples/performance.py5
-rw-r--r--examples/plot.py5
-rw-r--r--pyspike/DiscreteFunc.py64
-rw-r--r--pyspike/PieceWiseConstFunc.py13
-rw-r--r--pyspike/PieceWiseLinFunc.py17
-rw-r--r--pyspike/SpikeTrain.py15
-rw-r--r--pyspike/__init__.py25
-rw-r--r--pyspike/cython/cython_distances.pyx55
-rw-r--r--pyspike/cython/cython_profiles.pyx61
-rw-r--r--pyspike/cython/python_backend.py53
-rw-r--r--pyspike/generic.py18
-rw-r--r--pyspike/isi_distance.py14
-rw-r--r--pyspike/psth.py2
-rw-r--r--pyspike/spike_distance.py15
-rw-r--r--pyspike/spike_sync.py12
-rw-r--r--pyspike/spikes.py3
-rw-r--r--setup.py32
-rw-r--r--test/test_distance.py70
-rw-r--r--test/test_empty.py16
-rw-r--r--test/test_regression/test_regression_15.py80
-rw-r--r--test/test_spikes.py9
27 files changed, 434 insertions, 175 deletions
diff --git a/.travis.yml b/.travis.yml
index 1035775..1562d21 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,6 +2,9 @@ language: python
python:
- "2.6"
- "2.7"
+ - "3.3"
+ - "3.4"
+ - "3.5"
env:
- CYTHON_INSTALL="pip install -q cython"
- CYTHON_INSTALL=""
diff --git a/Changelog b/Changelog
index 519dd3b..2be5e52 100644
--- a/Changelog
+++ b/Changelog
@@ -1,3 +1,9 @@
+PySpike v0.4:
+ * Python 3 support (thanks to Igor Gnatenko)
+ * list interface to SpikeTrain class
+ * disable_backend_warning property
+ * several bugfixes
+
PySpike v0.3:
* addition of __version__ attribute
* restructured docs, Readme now only contains basic examples
diff --git a/Contributors.txt b/Contributors.txt
index 83563c7..512aa7f 100644
--- a/Contributors.txt
+++ b/Contributors.txt
@@ -1,5 +1,6 @@
Python/C Programming:
- Mario Mulansky
+- Igor Gnatenko
Scientific Methods:
- Thomas Kreuz
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000..aed0ae0
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,7 @@
+include *.rst
+include *.txt
+include pyspike/cython/*.c
+include directionality/cython/*.c
+recursive-include examples *.py *.txt
+recursive-include test *.py *.txt
+recursive-include doc *
diff --git a/doc/conf.py b/doc/conf.py
index 8011ea9..7e13eae 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -64,9 +64,9 @@ copyright = u'2014-2015, Mario Mulansky'
# built documents.
#
# The short X.Y version.
-version = '0.3'
+version = '0.4'
# The full version, including alpha/beta/rc tags.
-release = '0.3.0'
+release = '0.4.0'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
diff --git a/examples/multivariate.py b/examples/multivariate.py
index 53dbf0f..93f8516 100644
--- a/examples/multivariate.py
+++ b/examples/multivariate.py
@@ -23,8 +23,8 @@ spike_trains = spk.load_spike_trains_from_txt("PySpike_testdata.txt",
t_loading = time.clock()
print("Number of spike trains: %d" % len(spike_trains))
-num_of_spikes = sum([len(spike_trains[i].spikes)
- for i in xrange(len(spike_trains))])
+num_of_spikes = sum([len(spike_trains[i])
+ for i in range(len(spike_trains))])
print("Number of spikes: %d" % num_of_spikes)
# calculate the multivariate spike distance
diff --git a/examples/performance.py b/examples/performance.py
index 1c31e8f..ec6c830 100644
--- a/examples/performance.py
+++ b/examples/performance.py
@@ -14,6 +14,9 @@ from datetime import datetime
import cProfile
import pstats
+# in case you dont have the cython backends, disable the warnings as follows:
+# spk.disable_backend_warning = True
+
M = 100 # number of spike trains
r = 1.0 # rate of Poisson spike times
T = 1E3 # length of spike trains
@@ -23,7 +26,7 @@ print("%d spike trains with %d spikes" % (M, int(r*T)))
spike_trains = []
t_start = datetime.now()
-for i in xrange(M):
+for i in range(M):
spike_trains.append(spk.generate_poisson_spikes(r, T))
t_end = datetime.now()
runtime = (t_end-t_start).total_seconds()
diff --git a/examples/plot.py b/examples/plot.py
index 9670286..1922939 100644
--- a/examples/plot.py
+++ b/examples/plot.py
@@ -16,12 +16,13 @@ import matplotlib.pyplot as plt
import pyspike as spk
+
spike_trains = spk.load_spike_trains_from_txt("PySpike_testdata.txt",
edges=(0, 4000))
-# plot the spike time
+# plot the spike times
for (i, spike_train) in enumerate(spike_trains):
- plt.plot(spike_train.spikes, i*np.ones_like(spike_train.spikes), 'o')
+ plt.scatter(spike_train, i*np.ones_like(spike_train), marker='|')
f = spk.isi_profile(spike_trains[0], spike_trains[1])
x, y = f.get_plottable_data()
diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py
index 17153ee..fe97bc2 100644
--- a/pyspike/DiscreteFunc.py
+++ b/pyspike/DiscreteFunc.py
@@ -2,10 +2,11 @@
# Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
# Distributed under the BSD License
-from __future__ import print_function
+from __future__ import absolute_import, print_function
import numpy as np
import collections
+import pyspike
##############################################################
@@ -79,7 +80,7 @@ class DiscreteFunc(object):
expected_mp = (averaging_window_size+1) * int(self.mp[0])
y_plot = np.zeros_like(self.y)
# compute the values in a loop, could be done in cython if required
- for i in xrange(len(y_plot)):
+ for i in range(len(y_plot)):
if self.mp[i] >= expected_mp:
# the current value contains already all the wanted
@@ -136,9 +137,8 @@ class DiscreteFunc(object):
: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
+ value = 0.0
+ multiplicity = 0.0
def get_indices(ival):
""" Retuns the indeces surrounding the given interval"""
@@ -151,28 +151,32 @@ 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]))
-
- # check if interval is as sequence
- assert isinstance(interval, collections.Sequence), \
- "Invalid value for `interval`. None, Sequence or Tuple expected."
- # check if interval is a sequence of intervals
- 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]),
- np.sum(self.mp[start_ind:end_ind]))
+ value = 1.0 * np.sum(self.y[1:-1])
+ multiplicity = np.sum(self.mp[1:-1])
else:
- value = 0.0
- multiplicity = 0.0
- for ival in interval:
+ # check if interval is as sequence
+ assert isinstance(interval, collections.Sequence), \
+ "Invalid value for `interval`. None, Sequence or Tuple \
+expected."
+ # check if interval is a sequence of intervals
+ if not isinstance(interval[0], collections.Sequence):
# find the indices corresponding to the interval
- 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])
+ start_ind, end_ind = get_indices(interval)
+ value = np.sum(self.y[start_ind:end_ind])
+ multiplicity = np.sum(self.mp[start_ind:end_ind])
+ else:
+ for ival in interval:
+ # find the indices corresponding to the interval
+ 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])
+ if multiplicity == 0.0:
+ # empty profile, return spike sync of 1
+ value = 1.0
+ multiplicity = 1.0
return (value, multiplicity)
- def avrg(self, interval=None):
+ def avrg(self, interval=None, normalize=True):
""" Computes the average of the interval sequence:
:math:`a = 1/N \\sum f_n` where N is the number of intervals.
@@ -185,7 +189,10 @@ class DiscreteFunc(object):
:rtype: float
"""
val, mp = self.integral(interval)
- return val/mp
+ if normalize:
+ return val/mp
+ else:
+ return val
def add(self, f):
""" Adds another `DiscreteFunc` function to this function.
@@ -199,15 +206,16 @@ class DiscreteFunc(object):
# cython version
try:
- from cython.cython_add import add_discrete_function_cython as \
+ from .cython.cython_add import add_discrete_function_cython as \
add_discrete_function_impl
except ImportError:
- print("Warning: add_discrete_function_cython not found. Make \
+ if not(pyspike.disable_backend_warning):
+ print("Warning: add_discrete_function_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 add_discrete_function_python as \
+ from .cython.python_backend import add_discrete_function_python as \
add_discrete_function_impl
self.x, self.y, self.mp = \
@@ -236,7 +244,7 @@ def average_profile(profiles):
assert len(profiles) > 1
avrg_profile = profiles[0].copy()
- for i in xrange(1, len(profiles)):
+ for i in range(1, len(profiles)):
avrg_profile.add(profiles[i])
avrg_profile.mul_scalar(1.0/len(profiles)) # normalize
diff --git a/pyspike/PieceWiseConstFunc.py b/pyspike/PieceWiseConstFunc.py
index 2705443..5ce5f27 100644
--- a/pyspike/PieceWiseConstFunc.py
+++ b/pyspike/PieceWiseConstFunc.py
@@ -2,10 +2,11 @@
# Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
# Distributed under the BSD License
-from __future__ import print_function
+from __future__ import absolute_import, print_function
import numpy as np
import collections
+import pyspike
##############################################################
@@ -188,14 +189,16 @@ class PieceWiseConstFunc(object):
# cython version
try:
- from cython.cython_add import add_piece_wise_const_cython as \
+ from .cython.cython_add import add_piece_wise_const_cython as \
add_piece_wise_const_impl
except ImportError:
- print("Warning: add_piece_wise_const_cython not found. Make sure \
-that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \
+ if not(pyspike.disable_backend_warning):
+ print("Warning: add_piece_wise_const_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 add_piece_wise_const_python as \
+ from .cython.python_backend import add_piece_wise_const_python as \
add_piece_wise_const_impl
self.x, self.y = add_piece_wise_const_impl(self.x, self.y, f.x, f.y)
diff --git a/pyspike/PieceWiseLinFunc.py b/pyspike/PieceWiseLinFunc.py
index c0dd475..8145e63 100644
--- a/pyspike/PieceWiseLinFunc.py
+++ b/pyspike/PieceWiseLinFunc.py
@@ -2,10 +2,11 @@
# Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
# Distributed under the BSD License
-from __future__ import print_function
+from __future__ import absolute_import, print_function
import numpy as np
import collections
+import pyspike
##############################################################
@@ -221,20 +222,22 @@ class PieceWiseLinFunc:
assert self.x[-1] == f.x[-1], "The functions have different intervals"
# python implementation
- # from python_backend import add_piece_wise_lin_python
+ # from .python_backend import add_piece_wise_lin_python
# self.x, self.y1, self.y2 = add_piece_wise_lin_python(
# self.x, self.y1, self.y2, f.x, f.y1, f.y2)
# cython version
try:
- from cython.cython_add import add_piece_wise_lin_cython as \
+ from .cython.cython_add import add_piece_wise_lin_cython as \
add_piece_wise_lin_impl
except ImportError:
- print("Warning: add_piece_wise_lin_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.")
+ if not(pyspike.disable_backend_warning):
+ print("Warning: add_piece_wise_lin_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 add_piece_wise_lin_python as \
+ from .cython.python_backend import add_piece_wise_lin_python as \
add_piece_wise_lin_impl
self.x, self.y1, self.y2 = add_piece_wise_lin_impl(
diff --git a/pyspike/SpikeTrain.py b/pyspike/SpikeTrain.py
index 9127b60..4b59a5d 100644
--- a/pyspike/SpikeTrain.py
+++ b/pyspike/SpikeTrain.py
@@ -32,6 +32,21 @@ class SpikeTrain(object):
self.t_start = 0.0
self.t_end = float(edges)
+ def __getitem__(self, index):
+ """ Returns the time of the spike given by index.
+
+ :param index: Index of the spike.
+ :return: spike time.
+ """
+ return self.spikes[index]
+
+ def __len__(self):
+ """ Returns the number of spikes.
+
+ :return: Number of spikes.
+ """
+ return len(self.spikes)
+
def sort(self):
""" Sorts the spike times of this spike train using `np.sort`
"""
diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 3e836bd..069090b 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -4,27 +4,28 @@ Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
Distributed under the BSD License
"""
+from __future__ import absolute_import
+
__all__ = ["isi_distance", "spike_distance", "spike_sync", "psth",
"spikes", "SpikeTrain", "PieceWiseConstFunc", "PieceWiseLinFunc",
- "DiscreteFunc"]
+ "DiscreteFunc", "directionality"]
-from PieceWiseConstFunc import PieceWiseConstFunc
-from PieceWiseLinFunc import PieceWiseLinFunc
-from DiscreteFunc import DiscreteFunc
-from SpikeTrain import SpikeTrain
+from .PieceWiseConstFunc import PieceWiseConstFunc
+from .PieceWiseLinFunc import PieceWiseLinFunc
+from .DiscreteFunc import DiscreteFunc
+from .SpikeTrain import SpikeTrain
-from isi_distance import isi_profile, isi_distance, isi_profile_multi,\
+from .isi_distance import isi_profile, isi_distance, isi_profile_multi,\
isi_distance_multi, isi_distance_matrix
-from spike_distance import spike_profile, spike_distance, spike_profile_multi,\
+from .spike_distance import spike_profile, spike_distance, spike_profile_multi,\
spike_distance_multi, spike_distance_matrix
-from spike_sync import spike_sync_profile, spike_sync,\
+from .spike_sync import spike_sync_profile, spike_sync,\
spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix
-from psth import psth
+from .psth import psth
-from spikes import load_spike_trains_from_txt, spike_train_from_string, \
+from .spikes import load_spike_trains_from_txt, spike_train_from_string, \
merge_spike_trains, generate_poisson_spikes
-
# define the __version__ following
# http://stackoverflow.com/questions/17583443
from pkg_resources import get_distribution, DistributionNotFound
@@ -42,3 +43,5 @@ except DistributionNotFound:
__version__ = 'Please install this project with setup.py'
else:
__version__ = _dist.version
+
+disable_backend_warning = False
diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx
index c4f2349..41baa60 100644
--- a/pyspike/cython/cython_distances.pyx
+++ b/pyspike/cython/cython_distances.pyx
@@ -176,6 +176,8 @@ def spike_distance_cython(double[:] t1, double[:] t2,
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
+ cdef double[:] t_aux1 = np.empty(2)
+ cdef double[:] t_aux2 = np.empty(2)
spike_value = 0.0
@@ -184,19 +186,27 @@ def spike_distance_cython(double[:] t1, double[:] t2,
with nogil: # release the interpreter to allow multithreading
t_last = t_start
- t_p1 = t_start
- t_p2 = 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_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_start, t_end)
+ 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])
dt_p1 = dt_f1
- s1 = dt_p1*(t_f1-t_start)/isi1
+ # s1 = dt_p1*(t_f1-t_start)/isi1
+ s1 = dt_p1
index1 = -1
else:
t_f1 = t1[1]
- dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, 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]
s1 = dt_p1
@@ -204,14 +214,15 @@ def spike_distance_cython(double[:] t1, double[:] t2,
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_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])
- s2 = dt_p2*(t_f2-t_start)/isi2
+ # s2 = dt_p2*(t_f2-t_start)/isi2
+ s2 = dt_p2
index2 = -1
else:
t_f2 = t2[1]
- dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, 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]
s2 = dt_p2
@@ -233,7 +244,7 @@ def spike_distance_cython(double[:] t1, double[:] t2,
if index1 < N1-1:
t_f1 = t1[index1+1]
else:
- t_f1 = t_end
+ t_f1 = t_aux1[1]
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)
@@ -243,14 +254,16 @@ def spike_distance_cython(double[:] t1, double[:] t2,
# 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)
+ t_aux2[0], t_aux2[1])
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
+ # s1 = dt_p1*(t_end-t1[N1-1])/isi1
+ # Eero's correction: no adjustment
+ s1 = dt_p1
# 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):
@@ -264,7 +277,7 @@ def spike_distance_cython(double[:] t1, double[:] t2,
if index2 < N2-1:
t_f2 = t2[index2+1]
else:
- t_f2 = t_end
+ t_f2 = t_aux2[1]
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)
@@ -274,14 +287,16 @@ def spike_distance_cython(double[:] t1, double[:] t2,
# 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)
+ t_aux1[0], t_aux1[1])
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
+ # s2 = dt_p2*(t_end-t2[N2-1])/isi2
+ # Eero's correction: no adjustment
+ s2 = dt_p2
# 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
@@ -298,27 +313,27 @@ def spike_distance_cython(double[:] t1, double[:] t2,
if index1 < N1-1:
t_f1 = t1[index1+1]
dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
- t_start, t_end)
+ t_aux2[0], t_aux2[1])
isi1 = t_f1 - t_p1
else:
- t_f1 = t_end
+ t_f1 = t_aux1[1]
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)
+ t_aux1[0], t_aux1[1])
isi2 = t_f2 - t_p2
else:
- t_f2 = t_end
+ t_f2 = t_aux2[1]
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
+ 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
diff --git a/pyspike/cython/cython_profiles.pyx b/pyspike/cython/cython_profiles.pyx
index f9893eb..f08de6e 100644
--- a/pyspike/cython/cython_profiles.pyx
+++ b/pyspike/cython/cython_profiles.pyx
@@ -181,6 +181,8 @@ def spike_profile_cython(double[:] t1, double[:] t2,
cdef double[:] spike_events
cdef double[:] y_starts
cdef double[:] y_ends
+ cdef double[:] t_aux1 = np.empty(2)
+ cdef double[:] t_aux2 = np.empty(2)
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
@@ -189,6 +191,10 @@ 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
+
spike_events = np.empty(N1+N2+2)
y_starts = np.empty(len(spike_events)-1)
@@ -196,19 +202,27 @@ def spike_profile_cython(double[:] t1, double[:] t2,
with nogil: # release the interpreter to allow multithreading
spike_events[0] = t_start
- t_p1 = t_start
- t_p2 = 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_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_start, t_end)
+ 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])
dt_p1 = dt_f1
- s1 = dt_p1*(t_f1-t_start)/isi1
+ # s1 = dt_p1*(t_f1-t_start)/isi1
+ s1 = dt_p1
index1 = -1
else:
t_f1 = t1[1]
- dt_f1 = get_min_dist_cython(t_f1, t2, N2, 0, t_start, 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]
s1 = dt_p1
@@ -216,14 +230,15 @@ def spike_profile_cython(double[:] t1, double[:] t2,
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_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])
- s2 = dt_p2*(t_f2-t_start)/isi2
+ # s2 = dt_p2*(t_f2-t_start)/isi2
+ s2 = dt_p2
index2 = -1
else:
t_f2 = t2[1]
- dt_f2 = get_min_dist_cython(t_f2, t1, N1, 0, t_start, 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]
s2 = dt_p2
@@ -245,7 +260,7 @@ def spike_profile_cython(double[:] t1, double[:] t2,
if index1 < N1-1:
t_f1 = t1[index1+1]
else:
- t_f1 = t_end
+ t_f1 = t_aux1[1]
spike_events[index] = t_p1
s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
y_ends[index-1] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1,
@@ -253,14 +268,16 @@ def spike_profile_cython(double[:] t1, double[:] t2,
# 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)
+ t_aux2[0], t_aux2[1])
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
+ # s1 = dt_p1*(t_end-t1[N1-1])/isi1
+ # Eero's correction: no adjustment
+ s1 = dt_p1
# s2 is the same as above, thus we can compute y2 immediately
y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1,
isi2)
@@ -275,7 +292,7 @@ def spike_profile_cython(double[:] t1, double[:] t2,
if index2 < N2-1:
t_f2 = t2[index2+1]
else:
- t_f2 = t_end
+ t_f2 = t_aux2[1]
spike_events[index] = t_p2
s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
y_ends[index-1] = (s1*isi2 + s2*isi1) / isi_avrg_cython(isi1,
@@ -283,14 +300,16 @@ def spike_profile_cython(double[:] t1, double[:] t2,
# 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)
+ t_aux1[0], t_aux1[1])
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
+ # s2 = dt_p2*(t_end-t2[N2-1])/isi2
+ # Eero's correction: no adjustment
+ s2 = dt_p2
# s2 is the same as above, thus we can compute y2 immediately
y_starts[index] = (s1*isi2 + s2*isi1)/isi_avrg_cython(isi1, isi2)
else: # t_f1 == t_f2 - generate only one event
@@ -306,19 +325,19 @@ def spike_profile_cython(double[:] t1, double[:] t2,
if index1 < N1-1:
t_f1 = t1[index1+1]
dt_f1 = get_min_dist_cython(t_f1, t2, N2, index2,
- t_start, t_end)
+ t_aux2[0], t_aux2[1])
isi1 = t_f1 - t_p1
else:
- t_f1 = t_end
+ t_f1 = t_aux1[1]
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)
+ t_aux1[0], t_aux1[1])
isi2 = t_f2 - t_p2
else:
- t_f2 = t_end
+ t_f2 = t_aux2[1]
dt_f2 = dt_p2
isi2 = fmax(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
index += 1
@@ -330,8 +349,10 @@ def spike_profile_cython(double[:] t1, double[:] t2,
# 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*(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)
# end nogil
diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py
index 69a420f..a5f7ae4 100644
--- a/pyspike/cython/python_backend.py
+++ b/pyspike/cython/python_backend.py
@@ -144,35 +144,44 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
y_starts = np.empty(len(spike_events)-1)
y_ends = np.empty(len(spike_events)-1)
+ 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_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]
+
spike_events[0] = t_start
- t_p1 = t_start
- t_p2 = t_start
if t1[0] > t_start:
t_f1 = t1[0]
- dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end)
+ 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])
- s1 = dt_p1*(t_f1-t_start)/isi1
+ # s1 = dt_p1*(t_f1-t_start)/isi1
+ s1 = dt_p1
index1 = -1
else:
dt_p1 = 0.0
t_f1 = t1[1]
- dt_f1 = get_min_dist(t_f1, t2, 0, t_start, t_end)
+ dt_f1 = get_min_dist(t_f1, t2, 0, t_aux2[0], t_aux2[1])
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(t_f2, t1, 0, t_start, t_end)
+ 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])
- s2 = dt_p2*(t_f2-t_start)/isi2
+ # s2 = dt_p2*(t_f2-t_start)/isi2
+ s2 = dt_p2
index2 = -1
else:
dt_p2 = 0.0
t_f2 = t2[1]
- dt_f2 = get_min_dist(t_f2, t1, 0, t_start, t_end)
+ dt_f2 = get_min_dist(t_f2, t1, 0, t_aux1[0], t_aux1[1])
isi2 = t2[1]-t2[0]
s2 = dt_p2
index2 = 0
@@ -193,20 +202,22 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
if index1 < N1-1:
t_f1 = t1[index1+1]
else:
- t_f1 = t_end
+ t_f1 = t_aux1[1]
spike_events[index] = t_p1
s2 = (dt_p2*(t_f2-t_p1) + dt_f2*(t_p1-t_p2)) / isi2
y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
# now the next interval start value
if index1 < N1-1:
- dt_f1 = get_min_dist(t_f1, t2, index2, t_start, t_end)
+ dt_f1 = get_min_dist(t_f1, t2, index2, t_aux2[0], t_aux2[1])
isi1 = t_f1-t_p1
s1 = dt_p1
else:
dt_f1 = dt_p1
isi1 = max(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
+ # s1 = dt_p1*(t_end-t1[N1-1])/isi1
+ # Eero's correction: no adjustment
+ s1 = dt_p1
# s2 is the same as above, thus we can compute y2 immediately
y_starts[index] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
elif (index2 < N2-1) and (t_f1 > t_f2 or index1 == N1-1):
@@ -220,20 +231,22 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
if index2 < N2-1:
t_f2 = t2[index2+1]
else:
- t_f2 = t_end
+ t_f2 = t_aux2[1]
spike_events[index] = t_p2
s1 = (dt_p1*(t_f1-t_p2) + dt_f1*(t_p2-t_p1)) / isi1
y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
# now the next interval start value
if index2 < N2-1:
- dt_f2 = get_min_dist(t_f2, t1, index1, t_start, t_end)
+ dt_f2 = get_min_dist(t_f2, t1, index1, t_aux1[0], t_aux1[1])
isi2 = t_f2-t_p2
s2 = dt_p2
else:
dt_f2 = dt_p2
isi2 = max(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
+ # s2 = dt_p2*(t_end-t2[N2-1])/isi2
+ # Eero's adjustment: no correction
+ s2 = dt_p2
# s2 is the same as above, thus we can compute y2 immediately
y_starts[index] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)**2)
else: # t_f1 == t_f2 - generate only one event
@@ -248,18 +261,18 @@ def spike_distance_python(spikes1, spikes2, t_start, t_end):
y_starts[index] = 0.0
if index1 < N1-1:
t_f1 = t1[index1+1]
- dt_f1 = get_min_dist(t_f1, t2, index2, t_start, t_end)
+ dt_f1 = get_min_dist(t_f1, t2, index2, t_aux2[0], t_aux2[1])
isi1 = t_f1 - t_p1
else:
- t_f1 = t_end
+ t_f1 = t_aux1[1]
dt_f1 = dt_p1
isi1 = max(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(t_f2, t1, index1, t_start, t_end)
+ dt_f2 = get_min_dist(t_f2, t1, index1, t_aux1[0], t_aux1[1])
isi2 = t_f2 - t_p2
else:
- t_f2 = t_end
+ t_f2 = t_aux2[1]
dt_f2 = dt_p2
isi2 = max(t_end-t2[N2-1], t2[N2-1]-t2[N2-2])
index += 1
@@ -271,8 +284,8 @@ def spike_distance_python(spikes1, spikes2, t_start, 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 # *(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)
# use only the data added above
diff --git a/pyspike/generic.py b/pyspike/generic.py
index 41affcb..5ad06f1 100644
--- a/pyspike/generic.py
+++ b/pyspike/generic.py
@@ -7,6 +7,7 @@ Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
Distributed under the BSD License
"""
+from __future__ import division
import numpy as np
@@ -37,13 +38,15 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
"""
L1 = len(pairs1)
if L1 > 1:
- dist_prof1 = divide_and_conquer(pairs1[:L1/2], pairs1[L1/2:])
+ 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:])
+ 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]])
@@ -63,8 +66,8 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
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:])
+ 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]])
@@ -135,12 +138,13 @@ def _generic_distance_matrix(spike_trains, dist_function,
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:]]
+ pairs = [(i, j) for i in range(len(indices))
+ for j in range(i+1, len(indices))]
distance_matrix = np.zeros((len(indices), len(indices)))
for i, j in pairs:
- d = dist_function(spike_trains[i], spike_trains[j], interval)
+ d = dist_function(spike_trains[indices[i]], spike_trains[indices[j]],
+ interval)
distance_matrix[i, j] = d
distance_matrix[j, i] = d
return distance_matrix
diff --git a/pyspike/isi_distance.py b/pyspike/isi_distance.py
index 5ea555d..0ae7393 100644
--- a/pyspike/isi_distance.py
+++ b/pyspike/isi_distance.py
@@ -2,6 +2,9 @@
# Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
# Distributed under the BSD License
+from __future__ import absolute_import
+
+import pyspike
from pyspike import PieceWiseConstFunc
from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \
_generic_distance_matrix
@@ -31,14 +34,15 @@ def isi_profile(spike_train1, spike_train2):
# load cython implementation
try:
- from cython.cython_profiles import isi_profile_cython \
+ 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 \
+ if not(pyspike.disable_backend_warning):
+ print("Warning: isi_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 isi_distance_python \
+ from .cython.python_backend import isi_distance_python \
as isi_profile_impl
times, values = isi_profile_impl(spike_train1.get_spikes_non_empty(),
@@ -72,7 +76,7 @@ def isi_distance(spike_train1, spike_train2, interval=None):
# distance over the whole interval is requested: use specific function
# for optimal performance
try:
- from cython.cython_distances import isi_distance_cython \
+ from .cython.cython_distances import isi_distance_cython \
as isi_distance_impl
return isi_distance_impl(spike_train1.get_spikes_non_empty(),
diff --git a/pyspike/psth.py b/pyspike/psth.py
index 4027215..7cf1140 100644
--- a/pyspike/psth.py
+++ b/pyspike/psth.py
@@ -24,7 +24,7 @@ def psth(spike_trains, bin_size):
# N = len(spike_trains)
combined_spike_train = spike_trains[0].spikes
- for i in xrange(1, len(spike_trains)):
+ for i in range(1, len(spike_trains)):
combined_spike_train = np.append(combined_spike_train,
spike_trains[i].spikes)
diff --git a/pyspike/spike_distance.py b/pyspike/spike_distance.py
index ac2d260..e418283 100644
--- a/pyspike/spike_distance.py
+++ b/pyspike/spike_distance.py
@@ -2,6 +2,9 @@
# Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
# Distributed under the BSD License
+from __future__ import absolute_import
+
+import pyspike
from pyspike import PieceWiseLinFunc
from pyspike.generic import _generic_profile_multi, _generic_distance_multi, \
_generic_distance_matrix
@@ -31,14 +34,15 @@ def spike_profile(spike_train1, spike_train2):
# cython implementation
try:
- from cython.cython_profiles import spike_profile_cython \
+ from .cython.cython_profiles import spike_profile_cython \
as spike_profile_impl
except ImportError:
- print("Warning: spike_profile_cython not found. Make sure that \
+ if not(pyspike.disable_backend_warning):
+ 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 \
+ from .cython.python_backend import spike_distance_python \
as spike_profile_impl
times, y_starts, y_ends = spike_profile_impl(
@@ -54,7 +58,8 @@ Falling back to slow python backend.")
############################################################
def spike_distance(spike_train1, spike_train2, interval=None):
""" Computes the spike-distance :math:`D_S` of the given spike trains. The
- spike-distance is the integral over the isi distance profile :math:`S(t)`:
+ spike-distance is the integral over the spike distance profile
+ :math:`S(t)`:
.. math:: D_S = \int_{T_0}^{T_1} S(t) dt.
@@ -73,7 +78,7 @@ def spike_distance(spike_train1, spike_train2, interval=None):
# distance over the whole interval is requested: use specific function
# for optimal performance
try:
- from cython.cython_distances import spike_distance_cython \
+ 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(),
diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py
index 40d98d2..3dc29ff 100644
--- a/pyspike/spike_sync.py
+++ b/pyspike/spike_sync.py
@@ -3,8 +3,11 @@
# Copyright 2014-2015, Mario Mulansky <mario.mulansky@gmx.net>
# Distributed under the BSD License
+from __future__ import absolute_import
+
import numpy as np
from functools import partial
+import pyspike
from pyspike import DiscreteFunc
from pyspike.generic import _generic_profile_multi, _generic_distance_matrix
@@ -36,14 +39,15 @@ def spike_sync_profile(spike_train1, spike_train2, max_tau=None):
# cython implementation
try:
- from cython.cython_profiles import coincidence_profile_cython \
+ from .cython.cython_profiles import coincidence_profile_cython \
as coincidence_profile_impl
except ImportError:
- print("Warning: spike_distance_cython not found. Make sure that \
+ if not(pyspike.disable_backend_warning):
+ 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 \
+ from .cython.python_backend import coincidence_python \
as coincidence_profile_impl
if max_tau is None:
@@ -71,7 +75,7 @@ def _spike_sync_values(spike_train1, spike_train2, interval, max_tau):
# distance over the whole interval is requested: use specific function
# for optimal performance
try:
- from cython.cython_distances import coincidence_value_cython \
+ from .cython.cython_distances import coincidence_value_cython \
as coincidence_value_impl
if max_tau is None:
max_tau = 0.0
diff --git a/pyspike/spikes.py b/pyspike/spikes.py
index 35d8533..b18d7eb 100644
--- a/pyspike/spikes.py
+++ b/pyspike/spikes.py
@@ -28,7 +28,8 @@ def spike_train_from_string(s, edges, sep=' ', is_sorted=False):
# load_spike_trains_txt
############################################################
def load_spike_trains_from_txt(file_name, edges,
- separator=' ', comment='#', is_sorted=False):
+ separator=' ', comment='#', is_sorted=False,
+ ignore_empty_lines=True):
""" Loads a number of spike trains from a text file. Each line of the text
file should contain one spike train as a sequence of spike times separated
by `separator`. Empty lines as well as lines starting with `comment` are
diff --git a/setup.py b/setup.py
index d853cdf..8ef431a 100644
--- a/setup.py
+++ b/setup.py
@@ -33,23 +33,29 @@ ext_modules = []
if use_cython: # Cython is available, compile .pyx -> .c
ext_modules += [
- Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.pyx"]),
- Extension("pyspike.cython.cython_profiles", ["pyspike/cython/cython_profiles.pyx"]),
- Extension("pyspike.cython.cython_distances", ["pyspike/cython/cython_distances.pyx"]),
+ Extension("pyspike.cython.cython_add",
+ ["pyspike/cython/cython_add.pyx"]),
+ Extension("pyspike.cython.cython_profiles",
+ ["pyspike/cython/cython_profiles.pyx"]),
+ Extension("pyspike.cython.cython_distances",
+ ["pyspike/cython/cython_distances.pyx"])
]
cmdclass.update({'build_ext': build_ext})
elif use_c: # c files are there, compile to binaries
ext_modules += [
- Extension("pyspike.cython.cython_add", ["pyspike/cython/cython_add.c"]),
- Extension("pyspike.cython.cython_profiles", ["pyspike/cython/cython_profiles.c"]),
- Extension("pyspike.cython.cython_distances", ["pyspike/cython/cython_distances.c"]),
+ Extension("pyspike.cython.cython_add",
+ ["pyspike/cython/cython_add.c"]),
+ Extension("pyspike.cython.cython_profiles",
+ ["pyspike/cython/cython_profiles.c"]),
+ Extension("pyspike.cython.cython_distances",
+ ["pyspike/cython/cython_distances.c"])
]
# neither cython nor c files available -> automatic fall-back to python backend
setup(
name='pyspike',
packages=find_packages(exclude=['doc']),
- version='0.3.0',
+ version='0.4.0',
cmdclass=cmdclass,
ext_modules=ext_modules,
include_dirs=[numpy.get_include()],
@@ -78,10 +84,10 @@ train similarity',
'Programming Language :: Python :: 2',
'Programming Language :: Python :: 2.6',
'Programming Language :: Python :: 2.7',
- ],
- package_data={
- 'pyspike': ['cython/cython_add.c', 'cython/cython_profiles.c',
- 'cython_distances.c'],
- 'test': ['Spike_testdata.txt']
- }
+
+ 'Programming Language :: Python :: 3',
+ 'Programming Language :: Python :: 3.3',
+ 'Programming Language :: Python :: 3.4',
+ 'Programming Language :: Python :: 3.5',
+ ]
)
diff --git a/test/test_distance.py b/test/test_distance.py
index e45ac16..8cf81e2 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -17,6 +17,8 @@ from numpy.testing import assert_equal, assert_almost_equal, \
import pyspike as spk
from pyspike import SpikeTrain
+import os
+TEST_PATH = os.path.dirname(os.path.realpath(__file__))
def test_isi():
# generate two spike trains:
@@ -36,6 +38,7 @@ def test_isi():
f = spk.isi_profile(t1, t2)
# print("ISI: ", f.y)
+ print("ISI value:", expected_isi_val)
assert_equal(f.x, expected_times)
assert_array_almost_equal(f.y, expected_isi, decimal=15)
@@ -73,8 +76,19 @@ def test_spike():
assert_equal(f.x, expected_times)
- assert_almost_equal(f.avrg(), 1.6624149659863946e-01, decimal=15)
- assert_almost_equal(f.y2[-1], 0.1394558, decimal=6)
+ # from SPIKY:
+ y_all = np.array([0.000000000000000000, 0.555555555555555580,
+ 0.222222222222222210, 0.305555555555555580,
+ 0.255102040816326536, 0.000000000000000000,
+ 0.000000000000000000, 0.255102040816326536,
+ 0.255102040816326536, 0.285714285714285698,
+ 0.285714285714285698, 0.285714285714285698])
+
+ #assert_array_almost_equal(f.y1, y_all[::2])
+ assert_array_almost_equal(f.y2, y_all[1::2])
+
+ assert_almost_equal(f.avrg(), 0.186309523809523814, decimal=15)
+ assert_equal(spk.spike_distance(t1, t2), f.avrg())
t1 = SpikeTrain([0.2, 0.4, 0.6, 0.7], 1.0)
t2 = SpikeTrain([0.3, 0.45, 0.8, 0.9, 0.95], 1.0)
@@ -99,6 +113,8 @@ def test_spike():
(expected_y1+expected_y2)/2)
expected_spike_val /= (expected_times[-1]-expected_times[0])
+ print("SPIKE value:", expected_spike_val)
+
f = spk.spike_profile(t1, t2)
assert_equal(f.x, expected_times)
@@ -117,9 +133,14 @@ def test_spike():
# for left and right values
s1_r = np.array([0.1, (0.1*0.1+0.1*0.1)/0.2, 0.1, 0.0, 0.0, 0.0, 0.0])
s1_l = np.array([0.1, (0.1*0.1+0.1*0.1)/0.2, 0.1, 0.0, 0.0, 0.0, 0.0])
- s2_r = np.array([0.1*0.1/0.3, 0.1*0.3/0.3, 0.1*0.2/0.3,
+ # s2_r = np.array([0.1*0.1/0.3, 0.1*0.3/0.3, 0.1*0.2/0.3,
+ # 0.0, 0.1, 0.0, 0.0])
+ # s2_l = np.array([0.1*0.1/0.3, 0.1*0.1/0.3, 0.1*0.2/0.3, 0.0,
+ # 0.1, 0.0, 0.0])
+ # eero's edge correction:
+ s2_r = np.array([0.1, 0.1*0.3/0.3, 0.1*0.2/0.3,
0.0, 0.1, 0.0, 0.0])
- s2_l = np.array([0.1*0.1/0.3, 0.1*0.1/0.3, 0.1*0.2/0.3, 0.0,
+ s2_l = np.array([0.1, 0.1*0.3/0.3, 0.1*0.2/0.3, 0.0,
0.1, 0.0, 0.0])
isi1 = np.array([0.2, 0.2, 0.2, 0.2, 0.2, 0.4])
isi2 = np.array([0.3, 0.3, 0.3, 0.1, 0.1, 0.4])
@@ -275,8 +296,8 @@ def test_multi_spike_sync():
expected, decimal=15)
# multivariate regression test
- spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt",
- edges=[0, 4000])
+ spike_trains = spk.load_spike_trains_from_txt(
+ os.path.join(TEST_PATH, "SPIKE_Sync_Test.txt"), edges=[0, 4000])
# extract all spike times
spike_times = np.array([])
for st in spike_trains:
@@ -309,10 +330,10 @@ def check_dist_matrix(dist_func, dist_matrix_func):
f_matrix = dist_matrix_func(spike_trains)
# check zero diagonal
- for i in xrange(4):
+ for i in range(4):
assert_equal(0.0, f_matrix[i, i])
- for i in xrange(4):
- for j in xrange(i+1, 4):
+ for i in range(4):
+ for j in range(i+1, 4):
assert_equal(f_matrix[i, j], f_matrix[j, i])
assert_equal(f12, f_matrix[1, 0])
assert_equal(f13, f_matrix[2, 0])
@@ -345,15 +366,15 @@ def test_regression_spiky():
assert_equal(isi_profile.y, 0.1/1.1 * np.ones_like(isi_profile.y))
spike_dist = spk.spike_distance(st1, st2)
- assert_equal(spike_dist, 2.1105878248735391e-01)
+ assert_equal(spike_dist, 0.211058782487353908)
spike_sync = spk.spike_sync(st1, st2)
assert_equal(spike_sync, 8.6956521739130432e-01)
# multivariate check
- spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt",
- (0.0, 4000.0))
+ spike_trains = spk.load_spike_trains_from_txt(
+ os.path.join(TEST_PATH, "PySpike_testdata.txt"), (0.0, 4000.0))
isi_dist = spk.isi_distance_multi(spike_trains)
# get the full precision from SPIKY
assert_almost_equal(isi_dist, 0.17051816816999129656, decimal=15)
@@ -363,16 +384,35 @@ def test_regression_spiky():
spike_dist = spk.spike_distance_multi(spike_trains)
# get the full precision from SPIKY
- assert_almost_equal(spike_dist, 2.4432433330596512e-01, decimal=15)
+ assert_almost_equal(spike_dist, 0.25188056475463755, decimal=15)
spike_sync = spk.spike_sync_multi(spike_trains)
# get the full precision from SPIKY
assert_equal(spike_sync, 0.7183531505298066)
+ # Eero's edge correction example
+ st1 = SpikeTrain([0.5, 1.5, 2.5], 6.0)
+ st2 = SpikeTrain([3.5, 4.5, 5.5], 6.0)
+
+ f = spk.spike_profile(st1, st2)
+
+ expected_times = np.array([0.0, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.0])
+ y_all = np.array([0.271604938271605, 0.271604938271605, 0.271604938271605,
+ 0.617283950617284, 0.617283950617284, 0.444444444444444,
+ 0.285714285714286, 0.285714285714286, 0.444444444444444,
+ 0.617283950617284, 0.617283950617284, 0.271604938271605,
+ 0.271604938271605, 0.271604938271605])
+ expected_y1 = y_all[::2]
+ expected_y2 = y_all[1::2]
+
+ assert_equal(f.x, expected_times)
+ assert_array_almost_equal(f.y1, expected_y1, decimal=14)
+ assert_array_almost_equal(f.y2, expected_y2, decimal=14)
+
def test_multi_variate_subsets():
- spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt",
- (0.0, 4000.0))
+ spike_trains = spk.load_spike_trains_from_txt(
+ os.path.join(TEST_PATH, "PySpike_testdata.txt"), (0.0, 4000.0))
sub_set = [1, 3, 5, 7]
spike_trains_sub_set = [spike_trains[i] for i in sub_set]
diff --git a/test/test_empty.py b/test/test_empty.py
index 48be25d..5a0042f 100644
--- a/test/test_empty.py
+++ b/test/test_empty.py
@@ -70,8 +70,8 @@ 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)
- assert_almost_equal(d, 0.4*0.4*1.0/(0.4+1.0)**2 + 0.6*0.4*1.0/(0.6+1.0)**2,
- decimal=15)
+ d_expect = 0.4*0.4*1.0/(0.4+1.0)**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])
@@ -139,6 +139,18 @@ def test_spike_sync_empty():
assert_array_almost_equal(prof.x, [0.0, 0.2, 0.8, 1.0], decimal=15)
assert_array_almost_equal(prof.y, [0.0, 0.0, 0.0, 0.0], decimal=15)
+ # test with empty intervals
+ st1 = SpikeTrain([2.0, 5.0], [0, 10.0])
+ st2 = SpikeTrain([2.1, 7.0], [0, 10.0])
+ st3 = SpikeTrain([5.1, 6.0], [0, 10.0])
+ res = spk.spike_sync_profile(st1, st2).avrg(interval=[3.0, 4.0])
+ assert_equal(res, 1.0)
+ res = spk.spike_sync(st1, st2, interval=[3.0, 4.0])
+ assert_equal(res, 1.0)
+
+ sync_matrix = spk.spike_sync_matrix([st1, st2, st3], interval=[3.0, 4.0])
+ assert_array_equal(sync_matrix, np.ones((3, 3)) - np.diag(np.ones(3)))
+
if __name__ == "__main__":
test_get_non_empty()
diff --git a/test/test_regression/test_regression_15.py b/test/test_regression/test_regression_15.py
new file mode 100644
index 0000000..dcacae2
--- /dev/null
+++ b/test/test_regression/test_regression_15.py
@@ -0,0 +1,80 @@
+""" test_regression_15.py
+
+Regression test for Issue #15
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+
+from __future__ import division
+
+import numpy as np
+from numpy.testing import assert_equal, assert_almost_equal, \
+ assert_array_almost_equal
+
+import pyspike as spk
+
+import os
+TEST_PATH = os.path.dirname(os.path.realpath(__file__))
+TEST_DATA = os.path.join(TEST_PATH, "..", "SPIKE_Sync_Test.txt")
+
+def test_regression_15_isi():
+ # load spike trains
+ spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=[0, 4000])
+
+ N = len(spike_trains)
+
+ dist_mat = spk.isi_distance_matrix(spike_trains)
+ assert_equal(dist_mat.shape, (N, N))
+
+ ind = np.arange(N//2)
+ dist_mat = spk.isi_distance_matrix(spike_trains, ind)
+ assert_equal(dist_mat.shape, (N//2, N//2))
+
+ ind = np.arange(N//2, N)
+ dist_mat = spk.isi_distance_matrix(spike_trains, ind)
+ assert_equal(dist_mat.shape, (N//2, N//2))
+
+
+def test_regression_15_spike():
+ # load spike trains
+ spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=[0, 4000])
+
+ N = len(spike_trains)
+
+ dist_mat = spk.spike_distance_matrix(spike_trains)
+ assert_equal(dist_mat.shape, (N, N))
+
+ ind = np.arange(N//2)
+ dist_mat = spk.spike_distance_matrix(spike_trains, ind)
+ assert_equal(dist_mat.shape, (N//2, N//2))
+
+ ind = np.arange(N//2, N)
+ dist_mat = spk.spike_distance_matrix(spike_trains, ind)
+ assert_equal(dist_mat.shape, (N//2, N//2))
+
+
+def test_regression_15_sync():
+ # load spike trains
+ spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=[0, 4000])
+
+ N = len(spike_trains)
+
+ dist_mat = spk.spike_sync_matrix(spike_trains)
+ assert_equal(dist_mat.shape, (N, N))
+
+ ind = np.arange(N//2)
+ dist_mat = spk.spike_sync_matrix(spike_trains, ind)
+ assert_equal(dist_mat.shape, (N//2, N//2))
+
+ ind = np.arange(N//2, N)
+ dist_mat = spk.spike_sync_matrix(spike_trains, ind)
+ assert_equal(dist_mat.shape, (N//2, N//2))
+
+
+if __name__ == "__main__":
+ test_regression_15_isi()
+ test_regression_15_spike()
+ test_regression_15_sync()
diff --git a/test/test_spikes.py b/test/test_spikes.py
index d4eb131..609a819 100644
--- a/test/test_spikes.py
+++ b/test/test_spikes.py
@@ -13,10 +13,12 @@ from numpy.testing import assert_equal
import pyspike as spk
+import os
+TEST_PATH = os.path.dirname(os.path.realpath(__file__))
+TEST_DATA = os.path.join(TEST_PATH, "PySpike_testdata.txt")
def test_load_from_txt():
- spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt",
- edges=(0, 4000))
+ spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=(0, 4000))
assert len(spike_trains) == 40
# check the first spike train
@@ -48,8 +50,7 @@ def check_merged_spikes(merged_spikes, spike_trains):
def test_merge_spike_trains():
# first load the data
- spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt",
- edges=(0, 4000))
+ spike_trains = spk.load_spike_trains_from_txt(TEST_DATA, edges=(0, 4000))
merged_spikes = spk.merge_spike_trains([spike_trains[0], spike_trains[1]])
# test if result is sorted