summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Readme.rst94
-rw-r--r--examples/averages.py15
-rw-r--r--examples/distance_matrix.py5
-rw-r--r--examples/spike_sync.py39
-rw-r--r--pyspike/__init__.py7
-rw-r--r--pyspike/cython_add.pyx59
-rw-r--r--pyspike/cython_distance.pyx81
-rw-r--r--pyspike/distances.py146
-rw-r--r--pyspike/function.py165
-rw-r--r--pyspike/python_backend.py112
-rw-r--r--pyspike/spikes.py4
-rw-r--r--test/SPIKE_Sync_Test.txt100
-rw-r--r--test/test_distance.py67
-rw-r--r--test/test_function.py30
14 files changed, 814 insertions, 110 deletions
diff --git a/Readme.rst b/Readme.rst
index 4053995..bcf6fa9 100644
--- a/Readme.rst
+++ b/Readme.rst
@@ -5,8 +5,8 @@ PySpike
:target: https://travis-ci.org/mariomulansky/PySpike
PySpike is a Python library for the numerical analysis of spike train similarity.
-Its core functionality is the implementation of the bivariate ISI_ and SPIKE_ distance [#]_ [#]_.
-Additionally, it provides functions to compute multivariate SPIKE and ISI distances, as well as averaging and general spike train processing.
+Its core functionality is the implementation of the bivariate ISI_ and SPIKE_ distance [#]_ [#]_ as well as SPIKE-Synchronization_ [#]_.
+Additionally, it provides functions to compute multivariate profiles, distance matrices, as well as averaging and general spike train processing.
All computation intensive parts are implemented in C via cython_ to reach a competitive performance (factor 100-200 over plain Python).
PySpike provides the same fundamental functionality as the SPIKY_ framework for Matlab, which additionally contains spike-train generators, more spike train distance measures and many visualization routines.
@@ -17,6 +17,8 @@ All source codes are published under the BSD_License_.
.. [#] Kreuz T, Chicharro D, Houghton C, Andrzejak RG, Mormann F, *Monitoring spike train synchrony.* J Neurophysiol 109, 1457 (2013) `[pdf] <http://wwwold.fi.isc.cnr.it/users/thomas.kreuz/images/Kreuz_JNeurophysiol_2013_SPIKE-distance.pdf>`_
+.. [#] Kreuz T, Mulansky M and Bozanic N, *SPIKY: A graphical user interface for monitoring spike train synchrony*, tbp (2015)
+
Requirements and Installation
-----------------------------
@@ -162,32 +164,11 @@ If you are only interested in the scalar ISI-distance and not the profile, you c
where :code:`interval` is optional, as above, and if omitted the ISI-distance is computed for the complete spike trains.
-Furthermore, PySpike provides the :code:`average_profile` function that can be used to compute the average profile of a list of given :code:`PieceWiseConstFunc` instances.
-
-.. code:: python
-
- isi_profile1 = spk.isi_profile(spike_trains[0], spike_trains[1])
- isi_profile2 = spk.isi_profile(spike_trains[0], spike_trains[2])
- isi_profile3 = spk.isi_profile(spike_trains[1], spike_trains[2])
-
- avrg_profile = spk.average_profile([isi_profile1, isi_profile2, isi_profile3])
- x, y = avrg_profile.get_plottable_data()
- plt.plot(x, y, label="Average ISI profile")
-
-Note the difference between the :code:`average_profile` function, which returns a :code:`PieceWiseConstFunc` (or :code:`PieceWiseLinFunc`, see below), and the :code:`avrg` member function above, that computes the integral over the time profile resulting in a single value.
-So to obtain overall average ISI-distance of a list of ISI profiles you can first compute the average profile using :code:`average_profile` and the use
-
-.. code:: python
-
- avrg_isi = avrg_profile.avrg()
-
-to obtain the final, scalar average ISI distance of the whole set (see also "Computing multivariate distance" below).
-
SPIKE-distance
..............
-To compute for the spike distance you use the function :code:`spike_profile` instead of :code:`isi_profile` above.
+To compute for the spike distance profile you use the function :code:`spike_profile` instead of :code:`isi_profile` above.
But the general approach is very similar:
.. code:: python
@@ -217,22 +198,52 @@ Again, you can use
to compute the SPIKE distance directly, if you are not interested in the profile at all.
The parameter :code:`interval` is optional and if neglected the whole spike train is used.
-Furthmore, you can use the :code:`average_profile` function to compute an average profile of a list of SPIKE-profiles:
+
+
+SPIKE synchronization
+..............
+
+**Important note:**
+
+------------------------------
+
+ SPIKE-Synchronization measures *similarity*.
+ That means, a value of zero indicates absence of synchrony, while a value of one denotes the presence of synchrony.
+ This is exactly opposite to the other two measures: ISI- and SPIKE-distance.
+
+----------------------
+
+
+SPIKE synchronization is another approach to measure spike synchrony.
+In contrast to the SPIKE- and ISI-distance, it measures similarity instead of dissimilarity, i.e. higher values represent larger synchrony.
+Another difference is that the SPIKE synchronization profile is only defined exactly at the spike times, not for the whole interval of the spike trains.
+Therefore, it is represented by a :code:`DiscreteFunction`.
+
+To compute for the spike synchronization profile, PySpike provides the function :code:`spike_sync_profile`.
+The general handling of the profile, however, is similar to the other profiles above:
.. code:: python
+
+ import matplotlib.pyplot as plt
+ import pyspike as spk
- avrg_profile = spk.average_profile([spike_profile1, spike_profile2,
- spike_profile3])
- x, y = avrg_profile.get_plottable_data()
- plt.plot(x, y, label="Average SPIKE profile")
+ spike_trains = spk.load_spike_trains_from_txt("PySpike_testdata.txt",
+ time_interval=(0, 4000))
+ spike_profile = spk.spike_sync_profile(spike_trains[0], spike_trains[1])
+ x, y = spike_profile.get_plottable_data()
+
+For the direct computation of the overall spike synchronization value within some interval, the :code:`spike_sync` function can be used:
+
+.. code:: python
+
+ spike_sync = spk.spike_sync(spike_trains[0], spike_trains[1], interval)
Computing multivariate profiles and distances
----------------------------------------------
-To compute the multivariate ISI- or SPIKE-profile of a set of spike trains, you can compute all bivariate profiles separately and then use the :code:`average_profile` function above.
-However, PySpike provides convenience functions for that purpose.
-The following example computes the multivariate ISI- and SPIKE-profile for a list of spike trains:
+To compute the multivariate ISI-profile, SPIKE-profile or SPIKE-Synchronization profile f a set of spike trains, PySpike provides multi-variate version of the profile function.
+The following example computes the multivariate ISI-, SPIKE- and SPIKE-Sync-profile for a list of spike trains:
.. code:: python
@@ -240,15 +251,16 @@ The following example computes the multivariate ISI- and SPIKE-profile for a lis
time_interval=(0, 4000))
avrg_isi_profile = spk.isi_profile_multi(spike_trains)
avrg_spike_profile = spk.spike_profile_multi(spike_trains)
+ avrg_spike_sync_profile = spk.spike_sync_profile_multi(spike_trains)
-Both functions take an optional parameter :code:`indices`, a list of indices that allows to define the spike trains that should be used for the multivariate profile.
-As before, if you are only interested in the distance values, and not in the profile, PySpike offers the functions: :code:`isi_distance_multi` and :code:`spike_distance_multi`, that return the scalar multivariate ISI- and SPIKE-distance.
-Both distance functions also accept an :code:`interval` parameter that can be used to specify the begin and end of the averaging interval as a pair of floats, if neglected the complete interval is used.
+All functions take an optional parameter :code:`indices`, a list of indices that allows to define the spike trains that should be used for the multivariate profile.
+As before, if you are only interested in the distance values, and not in the profile, PySpike offers the functions: :code:`isi_distance_multi`, :code:`spike_distance_multi` and :code:`spike_sync_multi`, that return the scalar overall multivariate ISI- and SPIKE-distance as well as the SPIKE-Synchronization value.
+Those functions also accept an :code:`interval` parameter that can be used to specify the begin and end of the averaging interval as a pair of floats, if neglected the complete interval is used.
Another option to characterize large sets of spike trains are distance matrices.
-Each entry in the distance matrix represents a bivariate distance of two spike trains.
-The distance matrix is symmetric and has zero values at the diagonal.
-The following example computes and plots the ISI- and SPIKE-distance matrix, where for the latter one the averaging is performed only the time interval T=0..1000.
+Each entry in the distance matrix represents a bivariate distance (similarity for SPIKE-Synchronization) of two spike trains.
+The distance matrix is symmetric and has zero values (ones) at the diagonal.
+The following example computes and plots the ISI- and SPIKE-distance matrix as well as the SPIKE-Synchronization-matrix, with different intervals.
.. code:: python
@@ -264,6 +276,11 @@ The following example computes and plots the ISI- and SPIKE-distance matrix, whe
plt.imshow(spike_distance, interpolation='none')
plt.title("SPIKE-distance")
+ plt.figure()
+ spike_sync = spk.spike_sync_matrix(spike_trains, interval=(2000,4000))
+ plt.imshow(spike_sync, interpolation='none')
+ plt.title("SPIKE-Sync")
+
plt.show()
@@ -284,6 +301,7 @@ Curie Initial Training Network* `Neural Engineering Transformative Technologies
.. _ISI: http://www.scholarpedia.org/article/Measures_of_spike_train_synchrony#ISI-distance
.. _SPIKE: http://www.scholarpedia.org/article/SPIKE-distance
+.. _SPIKE-Synchronization: http://www.scholarpedia.org/article/Measures_of_spike_train_synchrony#SPIKE_synchronization
.. _cython: http://www.cython.org
.. _SPIKY: http://wwwold.fi.isc.cnr.it/users/thomas.kreuz/Source-Code/SPIKY.html
.. _BSD_License: http://opensource.org/licenses/BSD-2-Clause
diff --git a/examples/averages.py b/examples/averages.py
index 6413ab2..c3e81e2 100644
--- a/examples/averages.py
+++ b/examples/averages.py
@@ -42,3 +42,18 @@ print("SPIKE-distance (0-1000): %.8f" % spike1)
print("SPIKE-distance (1000-2000): %.8f" % spike2)
print("SPIKE-distance (0-1000) and (2000-3000): %.8f" % spike3)
print("SPIKE-distance (1000-2000) and (3000-4000): %.8f" % spike4)
+print()
+
+f = spk.spike_sync_profile(spike_trains[0], spike_trains[1])
+
+print("SPIKE-Synchronization: %.8f" % f.avrg())
+
+spike_sync1 = f.avrg(interval=(0, 1000))
+spike_sync2 = f.avrg(interval=(1000, 2000))
+spike_sync3 = f.avrg(interval=[(0, 1000), (2000, 3000)])
+spike_sync4 = f.avrg(interval=[(1000, 2000), (3000, 4000)])
+
+print("SPIKE-Sync (0-1000): %.8f" % spike_sync1)
+print("SPIKE-Sync (1000-2000): %.8f" % spike_sync2)
+print("SPIKE-Sync (0-1000) and (2000-3000): %.8f" % spike_sync3)
+print("SPIKE-Sync (1000-2000) and (3000-4000): %.8f" % spike_sync4)
diff --git a/examples/distance_matrix.py b/examples/distance_matrix.py
index 142db2c..0e4d825 100644
--- a/examples/distance_matrix.py
+++ b/examples/distance_matrix.py
@@ -30,4 +30,9 @@ spike_distance = spk.spike_distance_matrix(spike_trains, interval=(0, 1000))
plt.imshow(spike_distance, interpolation='none')
plt.title("SPIKE-distance, T=0-1000")
+plt.figure()
+spike_sync = spk.spike_sync_matrix(spike_trains, interval=(2000, 4000))
+plt.imshow(spike_sync, interpolation='none')
+plt.title("SPIKE-Sync, T=2000-4000")
+
plt.show()
diff --git a/examples/spike_sync.py b/examples/spike_sync.py
new file mode 100644
index 0000000..a72dbbf
--- /dev/null
+++ b/examples/spike_sync.py
@@ -0,0 +1,39 @@
+from __future__ import print_function
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import pyspike as spk
+
+spike_trains = spk.load_spike_trains_from_txt("../test/SPIKE_Sync_Test.txt",
+ time_interval=(0, 4000))
+
+plt.figure()
+
+f = spk.spike_sync_profile(spike_trains[0], spike_trains[1])
+x, y = f.get_plottable_data()
+plt.plot(x, y, '--ok', label="SPIKE-SYNC profile")
+print(f.x)
+print(f.y)
+print(f.mp)
+
+print("Average:", f.avrg())
+
+
+f = spk.spike_profile(spike_trains[0], spike_trains[1])
+x, y = f.get_plottable_data()
+
+plt.plot(x, y, '-b', label="SPIKE-profile")
+
+plt.axis([0, 4000, -0.1, 1.1])
+plt.legend(loc="center right")
+
+plt.figure()
+
+f = spk.spike_sync_profile_multi(spike_trains)
+x, y = f.get_plottable_data()
+plt.plot(x, y, '-k', label="SPIKE-Sync profile")
+
+print("Average:", f.avrg())
+
+plt.show()
diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index dca5722..55687e6 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -6,12 +6,13 @@ Distributed under the BSD License
__all__ = ["function", "distances", "spikes"]
-from function import PieceWiseConstFunc, PieceWiseLinFunc, average_profile
+from function import PieceWiseConstFunc, PieceWiseLinFunc, \
+ DiscreteFunction, average_profile
from distances import isi_profile, isi_distance, \
spike_profile, spike_distance, \
- spike_sync_profile, \
+ spike_sync_profile, spike_sync, \
isi_profile_multi, isi_distance_multi, isi_distance_matrix, \
spike_profile_multi, spike_distance_multi, spike_distance_matrix, \
- spike_sync_profile_multi
+ spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix
from spikes import add_auxiliary_spikes, load_spike_trains_from_txt, \
spike_train_from_string, merge_spike_trains, generate_poisson_spikes
diff --git a/pyspike/cython_add.pyx b/pyspike/cython_add.pyx
index bfbe208..817799e 100644
--- a/pyspike/cython_add.pyx
+++ b/pyspike/cython_add.pyx
@@ -172,3 +172,62 @@ def add_piece_wise_lin_cython(double[:] x1, double[:] y11, double[:] y12,
return (np.array(x_new[:index+2]),
np.array(y1_new[:index+1]),
np.array(y2_new[:index+1]))
+
+
+############################################################
+# add_discrete_function_cython
+############################################################
+def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1,
+ double[:] x2, double[:] y2, double[:] mp2):
+
+ cdef double[:] x_new = np.empty(len(x1) + len(x2))
+ cdef double[:] y_new = np.empty_like(x_new)
+ cdef double[:] mp_new = np.empty_like(x_new)
+ cdef int index1 = 0
+ cdef int index2 = 0
+ cdef int index = 0
+ cdef int N1 = len(y1)
+ cdef int N2 = len(y2)
+ x_new[0] = x1[0]
+ while (index1+1 < N1) and (index2+1 < N2):
+ if x1[index1+1] < x2[index2+1]:
+ index1 += 1
+ index += 1
+ x_new[index] = x1[index1]
+ y_new[index] = y1[index1]
+ mp_new[index] = mp1[index1]
+ elif x1[index1+1] > x2[index2+1]:
+ index2 += 1
+ index += 1
+ x_new[index] = x2[index2]
+ y_new[index] = y2[index2]
+ mp_new[index] = mp2[index2]
+ else: # x1[index1+1] == x2[index2+1]
+ index1 += 1
+ index2 += 1
+ index += 1
+ x_new[index] = x1[index1]
+ y_new[index] = y1[index1] + y2[index2]
+ mp_new[index] = mp1[index1] + mp2[index2]
+ # one array reached the end -> copy the contents of the other to the end
+ if index1+1 < len(y1):
+ x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:]
+ y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:]
+ mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:]
+ index += len(x1)-index1-1
+ elif index2+1 < len(y2):
+ x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:]
+ y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:]
+ mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:]
+ index += len(x2)-index2-1
+ # else: # both arrays reached the end simultaneously
+ # x_new[index+1] = x1[-1]
+ # y_new[index+1] = y1[-1] + y2[-1]
+ # mp_new[index+1] = mp1[-1] + mp2[-1]
+
+ y_new[0] = y_new[1]
+ mp_new[0] = mp_new[1]
+
+ # the last value is again the end of the interval
+ # only use the data that was actually filled
+ return x_new[:index+1], y_new[:index+1], mp_new[:index+1]
diff --git a/pyspike/cython_distance.pyx b/pyspike/cython_distance.pyx
index 779ff94..489aab9 100644
--- a/pyspike/cython_distance.pyx
+++ b/pyspike/cython_distance.pyx
@@ -33,6 +33,7 @@ cimport numpy as np
from libc.math cimport fabs
from libc.math cimport fmax
+from libc.math cimport fmin
DTYPE = np.float
ctypedef np.float_t DTYPE_t
@@ -229,3 +230,83 @@ def spike_distance_cython(double[:] t1,
# 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]
+
+
+
+############################################################
+# coincidence_python
+############################################################
+cdef inline double get_tau(double[:] spikes1, double[:] spikes2, int i, int j):
+ cdef double m = 1E100 # some huge number
+ cdef int N1 = len(spikes1)-2
+ cdef int N2 = len(spikes2)-2
+ if i < N1:
+ m = fmin(m, spikes1[i+1]-spikes1[i])
+ if j < N2:
+ m = fmin(m, spikes2[j+1]-spikes2[j])
+ if i > 1:
+ m = fmin(m, spikes1[i]-spikes1[i-1])
+ if j > 1:
+ m = fmin(m, spikes2[j]-spikes2[j-1])
+ return 0.5*m
+
+
+############################################################
+# coincidence_cython
+############################################################
+def coincidence_cython(double[:] spikes1, double[:] spikes2):
+
+ cdef int N1 = len(spikes1)
+ cdef int N2 = len(spikes2)
+ cdef int i = 0
+ cdef int j = 0
+ cdef int n = 0
+ cdef double[:] st = np.zeros(N1 + N2 - 2) # spike times
+ cdef double[:] c = np.zeros(N1 + N2 - 2) # coincidences
+ cdef double[:] mp = np.ones(N1 + N2 - 2) # multiplicity
+ cdef double tau
+ while n < N1 + N2 - 2:
+ if spikes1[i+1] < spikes2[j+1]:
+ i += 1
+ n += 1
+ tau = get_tau(spikes1, spikes2, i, j)
+ st[n] = spikes1[i]
+ if j > 0 and spikes1[i]-spikes2[j] < tau:
+ # coincidence between the current spike and the previous spike
+ # both get marked with 1
+ c[n] = 1
+ c[n-1] = 1
+ elif spikes1[i+1] > spikes2[j+1]:
+ j += 1
+ n += 1
+ tau = get_tau(spikes1, spikes2, i, j)
+ st[n] = spikes2[j]
+ if i > 0 and spikes2[j]-spikes1[i] < tau:
+ # coincidence between the current spike and the previous spike
+ # both get marked with 1
+ c[n] = 1
+ c[n-1] = 1
+ else: # spikes1[i+1] = spikes2[j+1]
+ # advance in both spike trains
+ j += 1
+ i += 1
+ if i == N1-1 or j == N2-1:
+ break
+ n += 1
+ # add only one event, but with coincidence 2 and multiplicity 2
+ st[n] = spikes1[i]
+ c[n] = 2
+ mp[n] = 2
+
+ st = st[:n+2]
+ c = c[:n+2]
+ mp = mp[:n+2]
+
+ st[0] = spikes1[0]
+ st[len(st)-1] = spikes1[len(spikes1)-1]
+ c[0] = c[1]
+ c[len(c)-1] = c[len(c)-2]
+ mp[0] = mp[1]
+ mp[len(mp)-1] = mp[len(mp)-2]
+
+ return st, c, mp
diff --git a/pyspike/distances.py b/pyspike/distances.py
index 04ea77b..5476b6f 100644
--- a/pyspike/distances.py
+++ b/pyspike/distances.py
@@ -11,7 +11,7 @@ import numpy as np
import threading
from functools import partial
-from pyspike import PieceWiseConstFunc, PieceWiseLinFunc
+from pyspike import PieceWiseConstFunc, PieceWiseLinFunc, DiscreteFunction
############################################################
@@ -132,39 +132,57 @@ def spike_distance(spikes1, spikes2, interval=None):
############################################################
# spike_sync_profile
############################################################
-def spike_sync_profile(spikes1, spikes2, k=3):
+def spike_sync_profile(spikes1, spikes2):
+ """ Computes the spike-synchronization profile S_sync(t) of the two given
+ spike trains. Returns the profile as a DiscreteFunction object. The S_sync
+ values are either 1 or 0, indicating the presence or absence of a
+ coincidence. The spike trains are expected to have auxiliary spikes at the
+ beginning and end of the interval. Use the function add_auxiliary_spikes to
+ add those spikes to the spike train.
- assert k > 0
+ :param spikes1: ordered array of spike times with auxiliary spikes.
+ :param spikes2: ordered array of spike times with auxiliary spikes.
+ :returns: The spike-distance profile :math:`S_{sync}(t)`.
+ :rtype: :class:`pyspike.function.DiscreteFunction`
+
+ """
# cython implementation
try:
from cython_distance import coincidence_cython \
as coincidence_impl
except ImportError:
-# print("Warning: spike_distance_cython not found. Make sure that \
-# PySpike is installed by running\n 'python setup.py build_ext --inplace'!\n \
-# Falling back to slow python backend.")
+ 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 python_backend import coincidence_python \
as coincidence_impl
- st, c = coincidence_impl(spikes1, spikes2)
+ times, coincidences, multiplicity = coincidence_impl(spikes1, spikes2)
- dc = np.zeros(len(c))
- for i in xrange(2*k):
- dc[k:-k] += c[i:-2*k+i]
+ return DiscreteFunction(times, coincidences, multiplicity)
- for n in xrange(0, k):
- for i in xrange(n+k):
- dc[n] += c[i]
- dc[-n-1] += c[-i-1]
- for i in xrange(k-n-1):
- dc[n] += c[i]
- dc[-n-1] += c[-i-1]
- dc *= 1.0/k
+############################################################
+# spike_sync
+############################################################
+def spike_sync(spikes1, spikes2, interval=None):
+ """ Computes the spike synchronization value SYNC of the given spike
+ trains. The spike synchronization value is the computed as the total number
+ of coincidences divided by the total number of spikes:
- return PieceWiseConstFunc(st, dc)
+ .. math:: SYNC = \sum_n C_n / N.
+
+ :param spikes1: ordered array of spike times with auxiliary spikes.
+ :param spikes2: ordered array of spike times with auxiliary spikes.
+ :param interval: averaging interval given as a pair of floats (T0, T1),
+ if None the average over the whole function is computed.
+ :type interval: Pair of floats or None.
+ :returns: The spike synchronization value.
+ :rtype: double
+ """
+ return spike_sync_profile(spikes1, spikes2).avrg(interval)
############################################################
@@ -201,8 +219,7 @@ def _generic_profile_multi(spike_trains, pair_distance_func, indices=None):
for (i, j) in pairs[1:]:
current_dist = pair_distance_func(spike_trains[i], spike_trains[j])
average_dist.add(current_dist) # add to the average
- average_dist.mul_scalar(1.0/len(pairs)) # normalize
- return average_dist
+ return average_dist, len(pairs)
############################################################
@@ -273,7 +290,10 @@ def isi_profile_multi(spike_trains, indices=None):
:returns: The averaged isi profile :math:`<S_{isi}>(t)`
:rtype: :class:`pyspike.function.PieceWiseConstFunc`
"""
- return _generic_profile_multi(spike_trains, isi_profile, indices)
+ average_dist, M = _generic_profile_multi(spike_trains, isi_profile,
+ indices)
+ average_dist.mul_scalar(1.0/M) # normalize
+ return average_dist
############################################################
@@ -314,39 +334,65 @@ def spike_profile_multi(spike_trains, indices=None):
:rtype: :class:`pyspike.function.PieceWiseLinFunc`
"""
- return _generic_profile_multi(spike_trains, spike_profile, indices)
+ average_dist, M = _generic_profile_multi(spike_trains, spike_profile,
+ indices)
+ average_dist.mul_scalar(1.0/M) # normalize
+ return average_dist
+
+
+############################################################
+# spike_distance_multi
+############################################################
+def spike_distance_multi(spike_trains, indices=None, interval=None):
+ """ Computes the multi-variate spike distance for a set of spike trains.
+ That is the time average of the multi-variate spike profile:
+ S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt
+ where the sum goes over all pairs <i,j>
+
+ :param spike_trains: list of spike trains
+ :param indices: list of indices defining which spike trains to use,
+ if None all given spike trains are used (default=None)
+ :type indices: list or None
+ :param interval: averaging interval given as a pair of floats, if None
+ the average over the whole function is computed.
+ :type interval: Pair of floats or None.
+ :returns: The averaged spike distance S.
+ :rtype: double
+ """
+ return spike_profile_multi(spike_trains, indices).avrg(interval)
############################################################
# spike_profile_multi
############################################################
-def spike_sync_profile_multi(spike_trains, indices=None, k=3):
+def spike_sync_profile_multi(spike_trains, indices=None):
""" Computes the multi-variate spike synchronization profile for a set of
- spike trains. That is the average spike-distance of all pairs of spike
- trains:
- :math:`S_ss(t) = 2/((N(N-1)) sum_{<i,j>} S_{ss}^{i, j}`,
- where the sum goes over all pairs <i,j>
+ spike trains. For each spike in the set of spike trains, the multi-variate
+ profile is defined as the number of coincidences divided by the number of
+ spike trains pairs involving the spike train of containing this spike,
+ which is the number of spike trains minus one (N-1).
:param spike_trains: list of spike trains
:param indices: list of indices defining which spike trains to use,
if None all given spike trains are used (default=None)
:type indices: list or None
- :returns: The averaged spike profile :math:`<S_{ss}>(t)`
- :rtype: :class:`pyspike.function.PieceWiseConstFunc`
+ :returns: The multi-variate spike sync profile :math:`<S_{sync}>(t)`
+ :rtype: :class:`pyspike.function.DiscreteFunction`
"""
- prof_func = partial(spike_sync_profile, k=k)
- return _generic_profile_multi(spike_trains, prof_func, indices)
+ prof_func = partial(spike_sync_profile)
+ average_dist, M = _generic_profile_multi(spike_trains, prof_func,
+ indices)
+ # average_dist.mul_scalar(1.0/M) # no normalization here!
+ return average_dist
############################################################
# spike_distance_multi
############################################################
-def spike_distance_multi(spike_trains, indices=None, interval=None):
- """ Computes the multi-variate spike distance for a set of spike trains.
- That is the time average of the multi-variate spike profile:
- S_{spike} = \int_0^T 2/((N(N-1)) sum_{<i,j>} S_{spike}^{i, j} dt
- where the sum goes over all pairs <i,j>
+def spike_sync_multi(spike_trains, indices=None, interval=None):
+ """ Computes the multi-variate spike synchronization value for a set of
+ spike trains.
:param spike_trains: list of spike trains
:param indices: list of indices defining which spike trains to use,
@@ -355,10 +401,10 @@ def spike_distance_multi(spike_trains, indices=None, interval=None):
:param interval: averaging interval given as a pair of floats, if None
the average over the whole function is computed.
:type interval: Pair of floats or None.
- :returns: The averaged spike distance S.
+ :returns: The multi-variate spike synchronization value SYNC.
:rtype: double
"""
- return spike_profile_multi(spike_trains, indices).avrg(interval)
+ return spike_sync_profile_multi(spike_trains, indices).avrg(interval)
############################################################
@@ -434,3 +480,25 @@ def spike_distance_matrix(spike_trains, indices=None, interval=None):
"""
return _generic_distance_matrix(spike_trains, spike_distance,
indices, interval)
+
+
+############################################################
+# spike_sync_matrix
+############################################################
+def spike_sync_matrix(spike_trains, indices=None, interval=None):
+ """ Computes the overall spike-synchronization value of all pairs of
+ spike-trains.
+
+ :param spike_trains: list of spike trains
+ :param indices: list of indices defining which spike trains to use,
+ if None all given spike trains are used (default=None)
+ :type indices: list or None
+ :param interval: averaging interval given as a pair of floats, if None
+ the average over the whole function is computed.
+ :type interval: Pair of floats or None.
+ :returns: 2D array with the pair wise time spike synchronization values
+ :math:`SYNC_{ij}`
+ :rtype: np.array
+ """
+ return _generic_distance_matrix(spike_trains, spike_sync,
+ indices, interval)
diff --git a/pyspike/function.py b/pyspike/function.py
index 662606c..e0dadf6 100644
--- a/pyspike/function.py
+++ b/pyspike/function.py
@@ -136,15 +136,6 @@ class PieceWiseConstFunc(object):
a /= int_length
return a
- def avrg_function_value(self):
- """ Computes the average function value of the piece-wise const
- function: :math:`a = 1/N sum_i f_i` where N is the number of intervals.
-
- :returns: the average a.
- :rtype: float
- """
- return sum(self.y)/(len(self.y))
-
def add(self, f):
""" Adds another PieceWiseConst function to this function.
Note: only functions defined on the same interval can be summed.
@@ -364,6 +355,162 @@ that PySpike is installed by running\n 'python setup.py build_ext --inplace'! \
self.y2 *= fac
+##############################################################
+# DiscreteFunction
+##############################################################
+class DiscreteFunction(object):
+ """ A class representing values defined on a discrete set of points.
+ """
+
+ def __init__(self, x, y, multiplicity):
+ """ Constructs the discrete function.
+
+ :param x: array of length N defining the points at which the values are
+ defined.
+ :param y: array of length N degining the values at the points x.
+ :param multiplicity: array of length N defining the multiplicity of the
+ values.
+ """
+ # convert parameters to arrays, also ensures copying
+ self.x = np.array(x)
+ self.y = np.array(y)
+ self.mp = np.array(multiplicity)
+
+ def copy(self):
+ """ Returns a copy of itself
+
+ :rtype: :class:`DiscreteFunction`
+ """
+ return DiscreteFunction(self.x, self.y, self.mp)
+
+ def almost_equal(self, other, decimal=14):
+ """ Checks if the function is equal to another function up to `decimal`
+ precision.
+
+ :param other: another :class:`DiscreteFunction`
+ :returns: True if the two functions are equal up to `decimal` decimals,
+ False otherwise
+ :rtype: bool
+ """
+ eps = 10.0**(-decimal)
+ return np.allclose(self.x, other.x, atol=eps, rtol=0.0) and \
+ np.allclose(self.y, other.y, atol=eps, rtol=0.0) and \
+ np.allclose(self.mp, other.mp, atol=eps, rtol=0.0)
+
+ def get_plottable_data(self, k=0):
+ """ Returns two arrays containing x- and y-coordinates for immeditate
+ plotting of the interval sequence.
+
+ :returns: (x_plot, y_plot) containing plottable data
+ :rtype: pair of np.array
+
+ Example::
+
+ x, y = f.get_plottable_data()
+ plt.plot(x, y, '-o', label="Discrete function")
+ """
+
+ if k > 0:
+ raise NotImplementedError()
+
+ return 1.0*self.x, 1.0*self.y/self.mp
+
+ def integral(self, interval=None):
+ """ Returns the integral over the given interval. For the discrete
+ function, this amounts to the sum over all values divided by the total
+ multiplicity.
+
+ :param interval: integration interval given as a pair of floats, or a
+ sequence of pairs in case of multiple intervals, if
+ None the integral over the whole function is computed.
+ :type interval: Pair, sequence of pairs, or None.
+ :returns: the integral
+ :rtype: float
+ """
+
+ def get_indices(ival):
+ """ Retuns the indeces surrounding the given interval"""
+ start_ind = np.searchsorted(self.x, ival[0], side='right')
+ end_ind = np.searchsorted(self.x, ival[1], side='left')
+ assert start_ind > 0 and end_ind < len(self.x), \
+ "Invalid averaging interval"
+ return start_ind, end_ind
+
+ 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]))
+ else:
+ value = 0.0
+ multiplicity = 0.0
+ 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])
+ return value/multiplicity
+
+ def avrg(self, interval=None):
+ """ Computes the average of the interval sequence:
+ :math:`a = 1/N sum f_n ` where N is the number of intervals.
+
+ :param interval: averaging interval given as a pair of floats, a
+ sequence of pairs for averaging multiple intervals, or
+ None, if None the average over the whole function is
+ computed.
+ :type interval: Pair, sequence of pairs, or None.
+ :returns: the average a.
+ :rtype: float
+ """
+ return self.integral(interval)
+
+ def add(self, f):
+ """ Adds another `DiscreteFunction` function to this function.
+ Note: only functions defined on the same interval can be summed.
+
+ :param f: :class:`DiscreteFunction` function to be added.
+ :rtype: None
+ """
+ assert self.x[0] == f.x[0], "The functions have different intervals"
+ assert self.x[-1] == f.x[-1], "The functions have different intervals"
+
+ # cython version
+ try:
+ from cython_add import add_discrete_function_cython as \
+ add_discrete_function_impl
+ except ImportError:
+ 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 python_backend import add_discrete_function_python as \
+ add_discrete_function_impl
+
+ self.x, self.y, self.mp = \
+ add_discrete_function_impl(self.x, self.y, self.mp,
+ f.x, f.y, f.mp)
+
+ def mul_scalar(self, fac):
+ """ Multiplies the function with a scalar value
+
+ :param fac: Value to multiply
+ :type fac: double
+ :rtype: None
+ """
+ self.y *= fac
+
+
def average_profile(profiles):
""" Computes the average profile from the given ISI- or SPIKE-profiles.
diff --git a/pyspike/python_backend.py b/pyspike/python_backend.py
index 7f8ea8c..481daf9 100644
--- a/pyspike/python_backend.py
+++ b/pyspike/python_backend.py
@@ -248,52 +248,69 @@ def cumulative_sync_python(spikes1, spikes2):
def coincidence_python(spikes1, spikes2):
def get_tau(spikes1, spikes2, i, j):
- return 0.5*min([spikes1[i]-spikes1[i-1], spikes1[i+1]-spikes1[i],
- spikes2[j]-spikes2[j-1], spikes2[j+1]-spikes2[j]])
+ m = 1E100 # some huge number
+ if i < len(spikes1)-2:
+ m = min(m, spikes1[i+1]-spikes1[i])
+ if j < len(spikes2)-2:
+ m = min(m, spikes2[j+1]-spikes2[j])
+ if i > 1:
+ m = min(m, spikes1[i]-spikes1[i-1])
+ if j > 1:
+ m = min(m, spikes2[j]-spikes2[j-1])
+ return 0.5*m
N1 = len(spikes1)
N2 = len(spikes2)
i = 0
j = 0
n = 0
- st = np.zeros(N1 + N2 - 2)
- c = np.zeros(N1 + N2 - 3)
- c[0] = 0
- st[0] = 0
- while n < N1 + N2:
+ st = np.zeros(N1 + N2 - 2) # spike times
+ c = np.zeros(N1 + N2 - 2) # coincidences
+ mp = np.ones(N1 + N2 - 2) # multiplicity
+ while n < N1 + N2 - 2:
if spikes1[i+1] < spikes2[j+1]:
i += 1
n += 1
tau = get_tau(spikes1, spikes2, i, j)
st[n] = spikes1[i]
- if spikes1[i]-spikes2[j] > tau:
- c[n] = 0
- else:
+ if j > 0 and spikes1[i]-spikes2[j] < tau:
+ # coincidence between the current spike and the previous spike
+ # both get marked with 1
c[n] = 1
+ c[n-1] = 1
elif spikes1[i+1] > spikes2[j+1]:
j += 1
n += 1
tau = get_tau(spikes1, spikes2, i, j)
st[n] = spikes2[j]
- if spikes2[j]-spikes1[i] > tau:
- c[n] = 0
- else:
+ if i > 0 and spikes2[j]-spikes1[i] < tau:
+ # coincidence between the current spike and the previous spike
+ # both get marked with 1
c[n] = 1
+ c[n-1] = 1
else: # spikes1[i+1] = spikes2[j+1]
+ # advance in both spike trains
j += 1
i += 1
if i == N1-1 or j == N2-1:
break
n += 1
+ # add only one event, but with coincidence 2 and multiplicity 2
st[n] = spikes1[i]
- c[n] = 0
- n += 1
- st[n] = spikes1[i]
- c[n] = 1
- c[0] = c[2]
+ c[n] = 2
+ mp[n] = 2
+
+ st = st[:n+2]
+ c = c[:n+2]
+ mp = mp[:n+2]
+
st[0] = spikes1[0]
st[-1] = spikes1[-1]
+ c[0] = c[1]
+ c[-1] = c[-2]
+ mp[0] = mp[1]
+ mp[-1] = mp[-2]
- return st, c
+ return st, c, mp
############################################################
@@ -409,3 +426,60 @@ def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22):
y2_new[index] = y12[-1]+y22[-1]
# only use the data that was actually filled
return x_new[:index+2], y1_new[:index+1], y2_new[:index+1]
+
+
+############################################################
+# add_discrete_function_python
+############################################################
+def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2):
+
+ x_new = np.empty(len(x1) + len(x2))
+ y_new = np.empty_like(x_new)
+ mp_new = np.empty_like(x_new)
+ x_new[0] = x1[0]
+ index1 = 0
+ index2 = 0
+ index = 0
+ while (index1+1 < len(y1)) and (index2+1 < len(y2)):
+ if x1[index1+1] < x2[index2+1]:
+ index1 += 1
+ index += 1
+ x_new[index] = x1[index1]
+ y_new[index] = y1[index1]
+ mp_new[index] = mp1[index1]
+ elif x1[index1+1] > x2[index2+1]:
+ index2 += 1
+ index += 1
+ x_new[index] = x2[index2]
+ y_new[index] = y2[index2]
+ mp_new[index] = mp2[index2]
+ else: # x1[index1+1] == x2[index2+1]
+ index1 += 1
+ index2 += 1
+ index += 1
+ x_new[index] = x1[index1]
+ y_new[index] = y1[index1] + y2[index2]
+ mp_new[index] = mp1[index1] + mp2[index2]
+ # one array reached the end -> copy the contents of the other to the end
+ if index1+1 < len(y1):
+ x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:]
+ y_new[index+1:index+1+len(x1)-index1-1] = y1[index1+1:]
+ mp_new[index+1:index+1+len(x1)-index1-1] = mp1[index1+1:]
+ index += len(x1)-index1-1
+ elif index2+1 < len(y2):
+ x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:]
+ y_new[index+1:index+1+len(x2)-index2-1] = y2[index2+1:]
+ mp_new[index+1:index+1+len(x2)-index2-1] = mp2[index2+1:]
+ index += len(x2)-index2-1
+ # else: # both arrays reached the end simultaneously
+ # x_new[index+1] = x1[-1]
+ # y_new[index+1] = y1[-1] + y2[-1]
+ # mp_new[index+1] = mp1[-1] + mp2[-1]
+
+ y_new[0] = y_new[1]
+ mp_new[0] = mp_new[1]
+
+ # the last value is again the end of the interval
+ # only use the data that was actually filled
+ return x_new[:index+1], y_new[:index+1], mp_new[:index+1]
+
diff --git a/pyspike/spikes.py b/pyspike/spikes.py
index aa25c48..6a3353e 100644
--- a/pyspike/spikes.py
+++ b/pyspike/spikes.py
@@ -67,7 +67,7 @@ def spike_train_from_string(s, sep=' ', is_sorted=False):
# load_spike_trains_txt
############################################################
def load_spike_trains_from_txt(file_name, time_interval=None,
- separator=' ', comment='#', sort=True):
+ separator=' ', comment='#', is_sorted=False):
""" 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
@@ -94,7 +94,7 @@ def load_spike_trains_from_txt(file_name, time_interval=None,
for line in spike_file:
if len(line) > 1 and not line.startswith(comment):
# use only the lines with actual data and not commented
- spike_train = spike_train_from_string(line, separator, sort)
+ spike_train = spike_train_from_string(line, separator, is_sorted)
if time_interval is not None: # add auxil. spikes if times given
spike_train = add_auxiliary_spikes(spike_train, time_interval)
spike_trains.append(spike_train)
diff --git a/test/SPIKE_Sync_Test.txt b/test/SPIKE_Sync_Test.txt
new file mode 100644
index 0000000..b97f777
--- /dev/null
+++ b/test/SPIKE_Sync_Test.txt
@@ -0,0 +1,100 @@
+61.000000 171.000000 268.000000 278.000000 326.000000 373.000000 400.000000 577.000000 793.000000 796.000000 798.000000 936.000000 994.000000 1026.000000 1083.000000 1097.000000 1187.000000 1228.000000 1400.000000 1522.000000 1554.000000 1579.000000 1661.000000 1895.000000 2040.000000 2082.000000 2264.000000 2502.000000 2689.000000 2922.000000 3093.000000 3276.000000 3495.000000 3693.000000 3900.000000
+
+195.000000 400.000000 518.000000 522.000000 565.000000 569.000000 572.000000 630.000000 802.000000 938.000000 1095.000000 1198.000000 1222.000000 1316.000000 1319.000000 1328.000000 1382.000000 1505.000000 1631.000000 1662.000000 1676.000000 1708.000000 1803.000000 1947.000000 1999.000000 2129.000000 2332.000000 2466.000000 2726.000000 2896.000000 3102.000000 3316.000000 3505.000000 3707.000000 3900.000000
+
+45.000000 111.000000 282.000000 319.000000 366.000000 400.000000 570.000000 633.000000 673.000000 750.000000 796.000000 1014.000000 1096.000000 1167.000000 1180.000000 1237.000000 1341.000000 1524.000000 1571.000000 1574.000000 1590.000000 1610.000000 1832.000000 1869.000000 1949.000000 1968.000000 2353.000000 2497.000000 2713.000000 2868.000000 3095.000000 3302.000000 3525.000000 3704.000000 3900.000000
+
+60.000000 135.000000 204.000000 260.000000 297.000000 361.000000 364.000000 400.000000 438.000000 631.000000 787.000000 794.000000 908.000000 927.000000 1205.000000 1251.000000 1315.000000 1463.000000 1546.000000 1548.000000 1569.000000 1604.000000 1705.000000 1733.000000 1994.000000 2146.000000 2306.000000 2554.000000 2691.000000 2905.000000 3090.000000 3296.000000 3508.000000 3702.000000 3900.000000
+
+159.000000 186.000000 308.000000 331.000000 400.000000 624.000000 758.000000 805.000000 811.000000 876.000000 1018.000000 1122.000000 1193.000000 1308.000000 1354.000000 1524.000000 1550.000000 1672.000000 1728.000000 1738.000000 1899.000000 1919.000000 1980.000000 1991.000000 2050.000000 2124.000000 2308.000000 2450.000000 2703.000000 2876.000000 3096.000000 3288.000000 3494.000000 3709.000000 3900.000000
+
+23.000000 134.000000 278.000000 400.000000 443.000000 555.000000 589.000000 597.000000 750.000000 754.000000 807.000000 1207.000000 1302.000000 1386.000000 1390.000000 1391.000000 1411.000000 1467.000000 1564.000000 1621.000000 1672.000000 1707.000000 1752.000000 1889.000000 1983.000000 2026.000000 2277.000000 2527.000000 2685.000000 2903.000000 3124.000000 3326.000000 3499.000000 3695.000000 3900.000000
+
+65.000000 120.000000 142.000000 200.000000 331.000000 400.000000 683.000000 701.000000 707.000000 788.000000 795.000000 880.000000 948.000000 972.000000 1197.000000 1282.000000 1354.000000 1451.000000 1543.000000 1549.000000 1555.000000 1624.000000 1723.000000 1766.000000 2023.000000 2071.000000 2285.000000 2478.000000 2766.000000 2883.000000 3100.000000 3292.000000 3502.000000 3700.000000 3900.000000
+
+83.000000 122.000000 214.000000 354.000000 400.000000 505.000000 614.000000 621.000000 697.000000 788.000000 846.000000 871.000000 878.000000 1174.000000 1204.000000 1215.000000 1317.000000 1353.000000 1365.000000 1453.000000 1463.000000 1540.000000 1832.000000 2016.000000 2023.000000 2290.000000 2449.000000 2708.000000 2881.000000 3090.000000 3329.000000 3489.000000 3705.000000 3900.000000
+
+56.000000 62.000000 143.000000 226.000000 259.000000 400.000000 439.000000 441.000000 569.000000 572.000000 639.000000 697.000000 808.000000 1162.000000 1178.000000 1250.000000 1360.000000 1427.000000 1598.000000 1667.000000 1671.000000 1780.000000 1865.000000 1902.000000 1972.000000 2092.000000 2318.000000 2548.000000 2741.000000 2888.000000 3096.000000 3304.000000 3518.000000 3705.000000 3900.000000
+
+110.000000 116.000000 215.000000 400.000000 462.000000 542.000000 602.000000 614.000000 619.000000 795.000000 1166.000000 1196.000000 1240.000000 1252.000000 1268.000000 1295.000000 1405.000000 1561.000000 1597.000000 1725.000000 1750.000000 1759.000000 1877.000000 1948.000000 2053.000000 2119.000000 2339.000000 2527.000000 2672.000000 2874.000000 3137.000000 3312.000000 3488.000000 3698.000000 3900.000000
+
+141.000000 159.000000 332.000000 333.000000 400.000000 756.000000 801.000000 843.000000 1161.000000 1208.000000 1225.000000 1246.000000 1418.000000 1448.000000 1501.000000 1559.000000 1578.000000 1684.000000 1751.000000 1797.000000 1815.000000 1818.000000 1948.000000 1975.000000 1989.000000 2110.000000 2360.000000 2453.000000 2704.000000 2906.000000 3106.000000 3286.000000 3491.000000 3697.000000 3900.000000
+
+61.000000 145.000000 151.000000 340.000000 400.000000 642.000000 741.000000 801.000000 901.000000 912.000000 939.000000 1072.000000 1180.000000 1216.000000 1271.000000 1336.000000 1344.000000 1584.000000 1608.000000 1617.000000 1648.000000 1695.000000 1789.000000 1835.000000 2053.000000 2089.000000 2223.000000 2531.000000 2688.000000 2901.000000 3114.000000 3268.000000 3496.000000 3703.000000 3900.000000
+
+313.000000 378.000000 400.000000 548.000000 657.000000 691.000000 715.000000 728.000000 802.000000 813.000000 1092.000000 1203.000000 1237.000000 1388.000000 1562.000000 1566.000000 1573.000000 1659.000000 1781.000000 1788.000000 1821.000000 1825.000000 1835.000000 1985.000000 1993.000000 2152.000000 2308.000000 2492.000000 2681.000000 2890.000000 3101.000000 3305.000000 3492.000000 3694.000000 3900.000000
+
+244.000000 311.000000 392.000000 400.000000 431.000000 441.000000 562.000000 577.000000 632.000000 791.000000 818.000000 875.000000 1020.000000 1059.000000 1134.000000 1164.000000 1201.000000 1238.000000 1273.000000 1387.000000 1562.000000 1609.000000 1831.000000 1949.000000 1961.000000 2088.000000 2329.000000 2509.000000 2691.000000 2902.000000 3096.000000 3279.000000 3506.000000 3704.000000 3900.000000
+
+11.000000 111.000000 159.000000 277.000000 334.000000 400.000000 480.000000 646.000000 804.000000 1122.000000 1129.000000 1178.000000 1198.000000 1233.000000 1359.000000 1374.000000 1411.000000 1476.000000 1477.000000 1571.000000 1582.000000 1622.000000 1706.000000 1867.000000 1988.000000 2094.000000 2233.000000 2512.000000 2671.000000 2931.000000 3111.000000 3292.000000 3488.000000 3691.000000 3900.000000
+
+57.000000 114.000000 328.000000 400.000000 442.000000 582.000000 662.000000 752.000000 766.000000 795.000000 1035.000000 1115.000000 1204.000000 1242.000000 1261.000000 1277.000000 1295.000000 1300.000000 1333.000000 1398.000000 1571.000000 1594.000000 1743.000000 1765.000000 2076.000000 2094.000000 2319.000000 2518.000000 2683.000000 2933.000000 3109.000000 3317.000000 3492.000000 3696.000000 3900.000000
+
+92.000000 102.000000 111.000000 190.000000 400.000000 446.000000 478.000000 630.000000 631.000000 805.000000 823.000000 918.000000 985.000000 1199.000000 1209.000000 1217.000000 1355.000000 1466.000000 1503.000000 1563.000000 1582.000000 1636.000000 1819.000000 1944.000000 1977.000000 2014.000000 2359.000000 2428.000000 2728.000000 2868.000000 3101.000000 3296.000000 3509.000000 3708.000000 3900.000000
+
+34.000000 66.000000 70.000000 113.000000 135.000000 238.000000 284.000000 400.000000 528.000000 766.000000 805.000000 921.000000 994.000000 1045.000000 1137.000000 1180.000000 1193.000000 1481.000000 1625.000000 1660.000000 1699.000000 1764.000000 1809.000000 1861.000000 1967.000000 2095.000000 2267.000000 2518.000000 2719.000000 2885.000000 3081.000000 3252.000000 3484.000000 3705.000000 3900.000000
+
+65.000000 90.000000 123.000000 199.000000 330.000000 400.000000 805.000000 1005.000000 1035.000000 1044.000000 1064.000000 1138.000000 1155.000000 1205.000000 1217.000000 1248.000000 1318.000000 1345.000000 1403.000000 1567.000000 1609.000000 1781.000000 1875.000000 1929.000000 2024.000000 2140.000000 2258.000000 2477.000000 2747.000000 2890.000000 3120.000000 3325.000000 3510.000000 3708.000000 3900.000000
+
+70.000000 221.000000 280.000000 400.000000 489.000000 786.000000 1016.000000 1027.000000 1029.000000 1145.000000 1186.000000 1195.000000 1256.000000 1304.000000 1314.000000 1476.000000 1618.000000 1657.000000 1730.000000 1748.000000 1802.000000 1812.000000 1832.000000 1947.000000 1999.000000 2027.000000 2288.000000 2532.000000 2679.000000 2919.000000 3077.000000 3316.000000 3516.000000 3705.000000 3900.000000
+
+153.000000 400.000000 474.000000 532.000000 568.000000 693.000000 738.000000 798.000000 806.000000 949.000000 1077.000000 1083.000000 1098.000000 1169.000000 1172.000000 1192.000000 1517.000000 1530.000000 1538.000000 1560.000000 1582.000000 1699.000000 1981.000000 1982.000000 2171.000000 2312.000000 2475.000000 2680.000000 2887.000000 3119.000000 3300.000000 3502.000000 3701.000000 3900.000000
+
+92.000000 152.000000 164.000000 400.000000 520.000000 619.000000 621.000000 647.000000 648.000000 808.000000 853.000000 865.000000 920.000000 949.000000 1148.000000 1225.000000 1231.000000 1348.000000 1375.000000 1635.000000 1646.000000 1686.000000 1711.000000 2004.000000 2079.000000 2347.000000 2501.000000 2709.000000 2930.000000 3061.000000 3319.000000 3494.000000 3690.000000 3900.000000
+
+74.000000 103.000000 247.000000 265.000000 400.000000 495.000000 501.000000 534.000000 552.000000 557.000000 601.000000 604.000000 792.000000 1003.000000 1138.000000 1195.000000 1252.000000 1325.000000 1336.000000 1425.000000 1646.000000 1657.000000 1795.000000 1990.000000 1992.000000 2062.000000 2300.000000 2509.000000 2690.000000 2913.000000 3066.000000 3276.000000 3460.000000 3700.000000 3900.000000
+
+45.000000 90.000000 156.000000 400.000000 468.000000 523.000000 577.000000 583.000000 708.000000 797.000000 815.000000 1052.000000 1063.000000 1189.000000 1215.000000 1218.000000 1266.000000 1288.000000 1299.000000 1512.000000 1519.000000 1584.000000 1769.000000 1791.000000 1964.000000 2082.000000 2348.000000 2530.000000 2703.000000 2893.000000 3031.000000 3290.000000 3504.000000 3702.000000 3900.000000
+
+140.000000 269.000000 400.000000 475.000000 492.000000 520.000000 569.000000 645.000000 727.000000 794.000000 819.000000 834.000000 957.000000 1122.000000 1210.000000 1374.000000 1471.000000 1485.000000 1515.000000 1574.000000 1668.000000 1732.000000 1743.000000 1917.000000 2041.000000 2104.000000 2294.000000 2453.000000 2662.000000 2894.000000 3128.000000 3301.000000 3489.000000 3705.000000 3900.000000
+
+28.000000 96.000000 112.000000 400.000000 426.000000 477.000000 584.000000 763.000000 804.000000 815.000000 1089.000000 1175.000000 1218.000000 1366.000000 1394.000000 1506.000000 1553.000000 1564.000000 1592.000000 1712.000000 1755.000000 1788.000000 1814.000000 1816.000000 1997.000000 2072.000000 2345.000000 2487.000000 2741.000000 2881.000000 3074.000000 3310.000000 3521.000000 3707.000000 3900.000000
+
+215.000000 286.000000 400.000000 461.000000 488.000000 489.000000 768.000000 796.000000 885.000000 919.000000 1188.000000 1253.000000 1432.000000 1476.000000 1521.000000 1524.000000 1566.000000 1590.000000 1684.000000 1714.000000 1733.000000 1776.000000 1816.000000 1943.000000 2016.000000 2031.000000 2308.000000 2488.000000 2642.000000 2832.000000 3120.000000 3293.000000 3507.000000 3702.000000 3900.000000
+
+77.000000 229.000000 302.000000 369.000000 400.000000 401.000000 404.000000 418.000000 804.000000 1026.000000 1110.000000 1179.000000 1187.000000 1227.000000 1456.000000 1458.000000 1476.000000 1629.000000 1630.000000 1640.000000 1697.000000 1734.000000 1785.000000 1919.000000 1956.000000 2057.000000 2324.000000 2416.000000 2656.000000 2889.000000 3126.000000 3323.000000 3491.000000 3696.000000 3900.000000
+
+244.000000 302.000000 400.000000 455.000000 533.000000 562.000000 673.000000 748.000000 791.000000 1120.000000 1136.000000 1191.000000 1235.000000 1238.000000 1296.000000 1336.000000 1447.000000 1466.000000 1551.000000 1594.000000 1691.000000 1744.000000 1897.000000 1959.000000 2060.000000 2109.000000 2230.000000 2564.000000 2717.000000 2900.000000 3089.000000 3320.000000 3491.000000 3712.000000 3900.000000
+
+3.000000 196.000000 199.000000 320.000000 339.000000 358.000000 400.000000 495.000000 690.000000 737.000000 760.000000 791.000000 849.000000 1027.000000 1194.000000 1220.000000 1242.000000 1313.000000 1354.000000 1435.000000 1523.000000 1621.000000 1775.000000 1788.000000 1999.000000 2074.000000 2245.000000 2478.000000 2750.000000 2893.000000 3113.000000 3302.000000 3485.000000 3690.000000 3900.000000
+
+206.000000 234.000000 261.000000 277.000000 341.000000 374.000000 400.000000 465.000000 613.000000 672.000000 745.000000 793.000000 799.000000 917.000000 954.000000 1144.000000 1180.000000 1283.000000 1484.000000 1574.000000 1575.000000 1795.000000 1965.000000 1984.000000 2086.000000 2093.000000 2312.000000 2501.000000 2738.000000 2879.000000 3084.000000 3270.000000 3483.000000 3701.000000 3900.000000
+
+154.000000 314.000000 400.000000 611.000000 615.000000 795.000000 823.000000 869.000000 908.000000 938.000000 960.000000 1024.000000 1049.000000 1068.000000 1185.000000 1420.000000 1441.000000 1496.000000 1610.000000 1709.000000 1712.000000 1740.000000 1885.000000 1917.000000 1992.000000 2079.000000 2224.000000 2508.000000 2713.000000 2861.000000 3096.000000 3300.000000 3509.000000 3696.000000 3900.000000
+
+26.000000 51.000000 83.000000 121.000000 343.000000 400.000000 625.000000 695.000000 697.000000 783.000000 803.000000 933.000000 1014.000000 1135.000000 1158.000000 1210.000000 1548.000000 1589.000000 1662.000000 1663.000000 1674.000000 1677.000000 1733.000000 1801.000000 1978.000000 2027.000000 2276.000000 2477.000000 2687.000000 2946.000000 3108.000000 3293.000000 3503.000000 3702.000000 3900.000000
+
+21.000000 39.000000 125.000000 198.000000 254.000000 400.000000 456.000000 510.000000 806.000000 881.000000 920.000000 1000.000000 1046.000000 1067.000000 1129.000000 1143.000000 1188.000000 1438.000000 1552.000000 1603.000000 1754.000000 1761.000000 1943.000000 1960.000000 1980.000000 2068.000000 2246.000000 2544.000000 2731.000000 2923.000000 3060.000000 3271.000000 3517.000000 3700.000000 3900.000000
+
+166.000000 237.000000 295.000000 300.000000 319.000000 369.000000 400.000000 407.000000 413.000000 428.000000 439.000000 804.000000 831.000000 899.000000 971.000000 1164.000000 1199.000000 1259.000000 1331.000000 1497.000000 1564.000000 1832.000000 1881.000000 1915.000000 1970.000000 2189.000000 2271.000000 2482.000000 2742.000000 2863.000000 3116.000000 3293.000000 3492.000000 3705.000000 3900.000000
+
+298.000000 323.000000 400.000000 423.000000 526.000000 662.000000 799.000000 821.000000 830.000000 933.000000 989.000000 1190.000000 1200.000000 1227.000000 1251.000000 1306.000000 1543.000000 1574.000000 1589.000000 1690.000000 1697.000000 1849.000000 1938.000000 1951.000000 2027.000000 2059.000000 2315.000000 2456.000000 2703.000000 2944.000000 3103.000000 3307.000000 3497.000000 3693.000000 3900.000000
+
+60.000000 172.000000 400.000000 413.000000 420.000000 600.000000 660.000000 690.000000 752.000000 789.000000 951.000000 1056.000000 1176.000000 1201.000000 1290.000000 1440.000000 1450.000000 1456.000000 1638.000000 1653.000000 1703.000000 1710.000000 1730.000000 1856.000000 2006.000000 2082.000000 2296.000000 2383.000000 2693.000000 2887.000000 3091.000000 3299.000000 3485.000000 3691.000000 3900.000000
+
+20.000000 127.000000 326.000000 369.000000 400.000000 521.000000 588.000000 595.000000 700.000000 798.000000 799.000000 858.000000 913.000000 1101.000000 1193.000000 1379.000000 1432.000000 1440.000000 1482.000000 1486.000000 1575.000000 1577.000000 1792.000000 1820.000000 1957.000000 2097.000000 2309.000000 2493.000000 2639.000000 2854.000000 3109.000000 3294.000000 3488.000000 3713.000000 3900.000000
+
+65.000000 119.000000 362.000000 400.000000 779.000000 803.000000 804.000000 897.000000 938.000000 984.000000 1147.000000 1207.000000 1266.000000 1319.000000 1373.000000 1579.000000 1596.000000 1626.000000 1644.000000 1650.000000 1725.000000 1776.000000 1851.000000 1965.000000 2023.000000 2116.000000 2331.000000 2552.000000 2727.000000 2855.000000 3081.000000 3268.000000 3521.000000 3698.000000 3900.000000
+
+4.000000 10.000000 50.000000 124.000000 151.000000 169.000000 314.000000 317.000000 400.000000 474.000000 549.000000 630.000000 704.000000 798.000000 1030.000000 1144.000000 1155.000000 1188.000000 1345.000000 1390.000000 1428.000000 1603.000000 1867.000000 1902.000000 1922.000000 1995.000000 2290.000000 2431.000000 2679.000000 2886.000000 3092.000000 3305.000000 3501.000000 3704.000000 3900.000000
+
+31.000000 37.000000 44.000000 211.000000 400.000000 445.000000 454.000000 602.000000 641.000000 760.000000 802.000000 850.000000 945.000000 1079.000000 1104.000000 1149.000000 1201.000000 1305.000000 1537.000000 1568.000000 1613.000000 1702.000000 1805.000000 1958.000000 1969.000000 2112.000000 2300.000000 2532.000000 2680.000000 2952.000000 3124.000000 3303.000000 3500.000000 3695.000000 3900.000000
+
+43.000000 259.000000 276.000000 342.000000 362.000000 375.000000 380.000000 400.000000 674.000000 800.000000 804.000000 809.000000 882.000000 947.000000 952.000000 1219.000000 1351.000000 1504.000000 1568.000000 1593.000000 1720.000000 1752.000000 1871.000000 1961.000000 2022.000000 2046.000000 2254.000000 2486.000000 2651.000000 2868.000000 3103.000000 3278.000000 3482.000000 3708.000000 3900.000000
+
+1.000000 219.000000 227.000000 235.000000 241.000000 400.000000 606.000000 618.000000 645.000000 738.000000 797.000000 943.000000 1217.000000 1343.000000 1424.000000 1448.000000 1578.000000 1661.000000 1706.000000 1765.000000 1903.000000 1915.000000 1975.000000 1987.000000 2084.000000 2324.000000 2490.000000 2671.000000 2865.000000 3063.000000 3331.000000 3505.000000 3702.000000 3900.000000
+
+103.000000 109.000000 356.000000 357.000000 400.000000 501.000000 714.000000 788.000000 793.000000 810.000000 859.000000 974.000000 1109.000000 1172.000000 1238.000000 1252.000000 1291.000000 1319.000000 1479.000000 1559.000000 1598.000000 1678.000000 1753.000000 1768.000000 1940.000000 2100.000000 2331.000000 2600.000000 2758.000000 2889.000000 3073.000000 3292.000000 3487.000000 3707.000000 3900.000000
+
+234.000000 362.000000 388.000000 399.000000 400.000000 407.000000 452.000000 483.000000 692.000000 721.000000 797.000000 809.000000 863.000000 1216.000000 1227.000000 1338.000000 1445.000000 1473.000000 1536.000000 1596.000000 1608.000000 1619.000000 1914.000000 1990.000000 2052.000000 2117.000000 2316.000000 2488.000000 2682.000000 2918.000000 3104.000000 3299.000000 3506.000000 3696.000000 3900.000000
+
+31.000000 91.000000 400.000000 422.000000 545.000000 587.000000 751.000000 794.000000 828.000000 962.000000 963.000000 1032.000000 1073.000000 1166.000000 1174.000000 1188.000000 1320.000000 1423.000000 1462.000000 1589.000000 1625.000000 1677.000000 1706.000000 1939.000000 2023.000000 2103.000000 2292.000000 2507.000000 2745.000000 2921.000000 3088.000000 3297.000000 3506.000000 3698.000000 3900.000000
+
+35.000000 92.000000 237.000000 296.000000 400.000000 515.000000 601.000000 613.000000 798.000000 852.000000 1201.000000 1248.000000 1257.000000 1286.000000 1429.000000 1616.000000 1633.000000 1656.000000 1778.000000 1819.000000 1838.000000 1864.000000 1903.000000 1918.000000 1991.000000 2106.000000 2315.000000 2455.000000 2690.000000 2891.000000 3084.000000 3280.000000 3488.000000 3698.000000 3900.000000
+
+20.000000 25.000000 172.000000 223.000000 274.000000 295.000000 368.000000 372.000000 400.000000 493.000000 717.000000 775.000000 795.000000 1015.000000 1200.000000 1319.000000 1444.000000 1559.000000 1592.000000 1694.000000 1743.000000 1757.000000 1841.000000 1859.000000 2043.000000 2075.000000 2336.000000 2461.000000 2764.000000 2905.000000 3099.000000 3293.000000 3494.000000 3697.000000 3900.000000
+
+76.000000 120.000000 130.000000 209.000000 396.000000 400.000000 559.000000 572.000000 671.000000 726.000000 803.000000 907.000000 1011.000000 1128.000000 1208.000000 1232.000000 1321.000000 1337.000000 1531.000000 1600.000000 1702.000000 1777.000000 1824.000000 1862.000000 1988.000000 1999.000000 2352.000000 2537.000000 2750.000000 2957.000000 3102.000000 3291.000000 3503.000000 3701.000000 3900.000000
+
+135.000000 250.000000 302.000000 310.000000 393.000000 400.000000 417.000000 684.000000 730.000000 804.000000 981.000000 982.000000 1081.000000 1197.000000 1240.000000 1313.000000 1409.000000 1431.000000 1473.000000 1498.000000 1561.000000 1615.000000 1782.000000 1925.000000 1979.000000 2149.000000 2293.000000 2490.000000 2676.000000 2932.000000 3117.000000 3315.000000 3502.000000 3690.000000 3900.000000
+
diff --git a/test/test_distance.py b/test/test_distance.py
index 7be0d9b..3b4329c 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -130,6 +130,30 @@ def test_spike():
decimal=16)
+def test_spike_sync():
+ spikes1 = np.array([1.0, 2.0, 3.0])
+ spikes2 = np.array([2.1])
+ spikes1 = spk.add_auxiliary_spikes(spikes1, 4.0)
+ spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0)
+ assert_almost_equal(spk.spike_sync(spikes1, spikes2),
+ 0.5, decimal=16)
+
+ spikes2 = np.array([3.1])
+ spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0)
+ assert_almost_equal(spk.spike_sync(spikes1, spikes2),
+ 0.5, decimal=16)
+
+ spikes2 = np.array([1.1])
+ spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0)
+ assert_almost_equal(spk.spike_sync(spikes1, spikes2),
+ 0.5, decimal=16)
+
+ spikes2 = np.array([0.9])
+ spikes2 = spk.add_auxiliary_spikes(spikes2, 4.0)
+ assert_almost_equal(spk.spike_sync(spikes1, spikes2),
+ 0.5, decimal=16)
+
+
def check_multi_profile(profile_func, profile_func_multi):
# generate spike trains:
t1 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6, 0.7]), 1.0)
@@ -148,6 +172,10 @@ def check_multi_profile(profile_func, profile_func_multi):
f_multi = profile_func_multi(spike_trains, [0, 1])
assert f_multi.almost_equal(f12, decimal=14)
+ f_multi1 = profile_func_multi(spike_trains, [1, 2, 3])
+ f_multi2 = profile_func_multi(spike_trains[1:])
+ assert f_multi1.almost_equal(f_multi2, decimal=14)
+
f = copy(f12)
f.add(f13)
f.add(f23)
@@ -172,6 +200,41 @@ def test_multi_spike():
check_multi_profile(spk.spike_profile, spk.spike_profile_multi)
+def test_multi_spike_sync():
+ # some basic multivariate check
+ spikes1 = np.array([100, 300, 400, 405, 410, 500, 700, 800,
+ 805, 810, 815, 900], dtype=float)
+ spikes2 = np.array([100, 200, 205, 210, 295, 350, 400, 510,
+ 600, 605, 700, 910], dtype=float)
+ spikes3 = np.array([100, 180, 198, 295, 412, 420, 510, 640,
+ 695, 795, 820, 920], dtype=float)
+ spikes1 = spk.add_auxiliary_spikes(spikes1, 1000)
+ spikes2 = spk.add_auxiliary_spikes(spikes2, 1000)
+ spikes3 = spk.add_auxiliary_spikes(spikes3, 1000)
+ assert_almost_equal(spk.spike_sync(spikes1, spikes2),
+ 0.5, decimal=15)
+ assert_almost_equal(spk.spike_sync(spikes1, spikes3),
+ 0.5, decimal=15)
+ assert_almost_equal(spk.spike_sync(spikes2, spikes3),
+ 0.5, decimal=15)
+
+ f = spk.spike_sync_profile_multi([spikes1, spikes2, spikes3])
+ # hands on definition of the average multivariate spike synchronization
+ # expected = (f1.integral() + f2.integral() + f3.integral()) / \
+ # (np.sum(f1.mp[1:-1])+np.sum(f2.mp[1:-1])+np.sum(f3.mp[1:-1]))
+ expected = 0.5
+ assert_almost_equal(f.avrg(), expected, decimal=15)
+ assert_almost_equal(spk.spike_sync_multi([spikes1, spikes2, spikes3]),
+ expected, decimal=15)
+
+ # multivariate regression test
+ spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt",
+ time_interval=(0, 4000))
+ f = spk.spike_sync_profile_multi(spike_trains)
+ assert_equal(np.sum(f.y[1:-1]), 39932)
+ assert_equal(np.sum(f.mp[1:-1]), 85554)
+
+
def check_dist_matrix(dist_func, dist_matrix_func):
# generate spike trains:
t1 = spk.add_auxiliary_spikes(np.array([0.2, 0.4, 0.6, 0.7]), 1.0)
@@ -210,6 +273,10 @@ def test_spike_matrix():
check_dist_matrix(spk.spike_distance, spk.spike_distance_matrix)
+def test_spike_sync_matrix():
+ check_dist_matrix(spk.spike_sync, spk.spike_sync_matrix)
+
+
def test_regression_spiky():
spike_trains = spk.load_spike_trains_from_txt("test/PySpike_testdata.txt",
(0.0, 4000.0))
diff --git a/test/test_function.py b/test/test_function.py
index ba87db8..933fd2e 100644
--- a/test/test_function.py
+++ b/test/test_function.py
@@ -203,6 +203,36 @@ def test_pwl_avrg():
assert_array_almost_equal(f_avrg.y2, y2_expected, decimal=16)
+def test_df():
+ # testing discrete function
+ x = [0.0, 1.0, 2.0, 2.5, 4.0]
+ y = [0.0, 1.0, 1.0, 0.0, 1.0]
+ mp = [1.0, 2.0, 1.0, 2.0, 1.0]
+ f = spk.DiscreteFunction(x, y, mp)
+ xp, yp = f.get_plottable_data()
+
+ xp_expected = [0.0, 1.0, 2.0, 2.5, 4.0]
+ yp_expected = [0.0, 0.5, 1.0, 0.0, 1.0]
+ assert_array_almost_equal(xp, xp_expected, decimal=16)
+ assert_array_almost_equal(yp, yp_expected, decimal=16)
+
+ assert_almost_equal(f.avrg(), 2.0/5.0, decimal=16)
+
+ # interval averaging
+ a = f.avrg([0.5, 2.4])
+ assert_almost_equal(a, 2.0/3.0, decimal=16)
+ a = f.avrg([1.5, 3.5])
+ assert_almost_equal(a, 1.0/3.0, decimal=16)
+ a = f.avrg((0.9, 3.5))
+ assert_almost_equal(a, 2.0/5.0, decimal=16)
+ a = f.avrg([1.1, 4.0])
+ assert_almost_equal(a, 1.0/3.0, decimal=16)
+
+ # averaging over multiple intervals
+ a = f.avrg([(0.5, 1.5), (1.5, 2.6)])
+ assert_almost_equal(a, 2.0/5.0, decimal=16)
+
+
if __name__ == "__main__":
test_pwc()
test_pwc_add()
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-15 16:23:30 +0000