diff options
-rw-r--r-- | pyspike/DiscreteFunc.py | 9 | ||||
-rw-r--r-- | pyspike/__init__.py | 4 | ||||
-rw-r--r-- | pyspike/cython/cython_add.pyx | 33 | ||||
-rw-r--r-- | pyspike/cython/cython_distances.pyx | 5 | ||||
-rw-r--r-- | pyspike/cython/python_backend.py | 33 | ||||
-rw-r--r-- | pyspike/spike_sync.py | 10 | ||||
-rw-r--r-- | pyspike/spikes.py | 41 | ||||
-rw-r--r-- | test/numeric/test_regression_random_spikes.py | 27 | ||||
-rw-r--r-- | test/test_distance.py | 33 | ||||
-rw-r--r-- | test/test_save_load.py | 42 |
10 files changed, 176 insertions, 61 deletions
diff --git a/pyspike/DiscreteFunc.py b/pyspike/DiscreteFunc.py index fe97bc2..caad290 100644 --- a/pyspike/DiscreteFunc.py +++ b/pyspike/DiscreteFunc.py @@ -170,10 +170,6 @@ expected." 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, normalize=True): @@ -190,7 +186,10 @@ expected." """ val, mp = self.integral(interval) if normalize: - return val/mp + if mp > 0: + return val/mp + else: + return 1.0 else: return val diff --git a/pyspike/__init__.py b/pyspike/__init__.py index 1e879c4..08253fb 100644 --- a/pyspike/__init__.py +++ b/pyspike/__init__.py @@ -23,8 +23,8 @@ from .spike_sync import spike_sync_profile, spike_sync,\ spike_sync_profile_multi, spike_sync_multi, spike_sync_matrix from .psth import psth -from .spikes import load_spike_trains_from_txt, spike_train_from_string, \ - import_spike_trains_from_time_series, \ +from .spikes import load_spike_trains_from_txt, save_spike_trains_to_txt, \ + spike_train_from_string, import_spike_trains_from_time_series, \ merge_spike_trains, generate_poisson_spikes # define the __version__ following diff --git a/pyspike/cython/cython_add.pyx b/pyspike/cython/cython_add.pyx index 8da1e53..25f1181 100644 --- a/pyspike/cython/cython_add.pyx +++ b/pyspike/cython/cython_add.pyx @@ -182,8 +182,8 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1, cdef int index1 = 0 cdef int index2 = 0 cdef int index = 0 - cdef int N1 = len(y1) - cdef int N2 = len(y2) + cdef int N1 = len(y1)-1 + cdef int N2 = len(y2)-1 x_new[0] = x1[0] while (index1+1 < N1) and (index2+1 < N2): if x1[index1+1] < x2[index2+1]: @@ -206,20 +206,21 @@ def add_discrete_function_cython(double[:] x1, double[:] y1, double[:] mp1, 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] + if index1+1 < N1: + x_new[index+1:index+1+N1-index1] = x1[index1+1:] + y_new[index+1:index+1+N1-index1] = y1[index1+1:] + mp_new[index+1:index+1+N1-index1] = mp1[index1+1:] + index += N1-index1 + elif index2+1 < N2: + x_new[index+1:index+1+N2-index2] = x2[index2+1:] + y_new[index+1:index+1+N2-index2] = y2[index2+1:] + mp_new[index+1:index+1+N2-index2] = mp2[index2+1:] + index += N2-index2 + else: # both arrays reached the end simultaneously + x_new[index+1] = x1[index1+1] + y_new[index+1] = y1[index1+1] + y2[index2+1] + mp_new[index+1] = mp1[index1+1] + mp2[index2+1] + index += 1 y_new[0] = y_new[1] mp_new[0] = mp_new[1] diff --git a/pyspike/cython/cython_distances.pyx b/pyspike/cython/cython_distances.pyx index 7dc9cb9..ac5f226 100644 --- a/pyspike/cython/cython_distances.pyx +++ b/pyspike/cython/cython_distances.pyx @@ -428,9 +428,4 @@ def coincidence_value_cython(double[:] spikes1, double[:] spikes2, mp += 2 coinc += 2 - if coinc == 0 and mp == 0: - # empty spike trains -> spike sync = 1 by definition - coinc = 1 - mp = 1 - return coinc, mp diff --git a/pyspike/cython/python_backend.py b/pyspike/cython/python_backend.py index a1e3a65..6b7209a 100644 --- a/pyspike/cython/python_backend.py +++ b/pyspike/cython/python_backend.py @@ -560,7 +560,9 @@ def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2): index1 = 0 index2 = 0 index = 0 - while (index1+1 < len(y1)) and (index2+1 < len(y2)): + N1 = len(x1)-1 + N2 = len(x2)-1 + while (index1+1 < N1) and (index2+1 < N2): if x1[index1+1] < x2[index2+1]: index1 += 1 index += 1 @@ -581,20 +583,21 @@ def add_discrete_function_python(x1, y1, mp1, x2, y2, mp2): 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] + if index1+1 < N1: + x_new[index+1:index+1+N1-index1] = x1[index1+1:] + y_new[index+1:index+1+N1-index1] = y1[index1+1:] + mp_new[index+1:index+1+N1-index1] = mp1[index1+1:] + index += N1-index1 + elif index2+1 < N2: + x_new[index+1:index+1+N2-index2] = x2[index2+1:] + y_new[index+1:index+1+N2-index2] = y2[index2+1:] + mp_new[index+1:index+1+N2-index2] = mp2[index2+1:] + index += N2-index2 + 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] + index += 1 y_new[0] = y_new[1] mp_new[0] = mp_new[1] diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py index 617dd86..80f7805 100644 --- a/pyspike/spike_sync.py +++ b/pyspike/spike_sync.py @@ -215,7 +215,10 @@ def spike_sync_bi(spike_train1, spike_train2, interval=None, max_tau=None): """ c, mp = _spike_sync_values(spike_train1, spike_train2, interval, max_tau) - return 1.0*c/mp + if mp == 0: + return 1.0 + else: + return 1.0*c/mp ############################################################ @@ -257,7 +260,10 @@ def spike_sync_multi(spike_trains, indices=None, interval=None, max_tau=None): coincidence += c mp += m - return coincidence/mp + if mp == 0.0: + return 1.0 + else: + return coincidence/mp ############################################################ diff --git a/pyspike/spikes.py b/pyspike/spikes.py index 1bf474c..486a4a0 100644 --- a/pyspike/spikes.py +++ b/pyspike/spikes.py @@ -2,6 +2,7 @@ # Copyright 2014, Mario Mulansky <mario.mulansky@gmx.net> # Distributed under the BSD License + import numpy as np from pyspike import SpikeTrain @@ -25,7 +26,7 @@ def spike_train_from_string(s, edges, sep=' ', is_sorted=False): ############################################################ -# load_spike_trains_txt +# load_spike_trains_from_txt ############################################################ def load_spike_trains_from_txt(file_name, edges, separator=' ', comment='#', is_sorted=False, @@ -47,13 +48,17 @@ def load_spike_trains_from_txt(file_name, edges, :returns: list of :class:`.SpikeTrain` """ spike_trains = [] - spike_file = open(file_name, 'r') - 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, edges, - separator, is_sorted) - spike_trains.append(spike_train) + with open(file_name, 'r') as spike_file: + for line in spike_file: + if not line.startswith(comment): # ignore comments + if len(line) > 1: + # ignore empty lines + spike_train = spike_train_from_string(line, edges, + separator, is_sorted) + spike_trains.append(spike_train) + elif not(ignore_empty_lines): + # add empty spike train + spike_trains.append(SpikeTrain([], edges)) return spike_trains @@ -83,6 +88,26 @@ def import_spike_trains_from_time_series(file_name, start_time, time_bin, ############################################################ +# save_spike_trains_to_txt +############################################################ +def save_spike_trains_to_txt(spike_trains, file_name, + separator=' ', precision=8): + """ Saves the given spike trains into a file with the given file name. + Each spike train will be stored in one line in the text file with the times + separated by `separator`. + + :param spike_trains: List of :class:`.SpikeTrain` objects + :param file_name: The name of the text file. + """ + # format string to print the spike times with given precision + format_str = "{0:.%de}" % precision + with open(file_name, 'w') as spike_file: + for st in spike_trains: + s = separator.join(map(format_str.format, st.spikes)) + spike_file.write(s+'\n') + + +############################################################ # merge_spike_trains ############################################################ def merge_spike_trains(spike_trains): diff --git a/test/numeric/test_regression_random_spikes.py b/test/numeric/test_regression_random_spikes.py index 6156bb4..73d39a6 100644 --- a/test/numeric/test_regression_random_spikes.py +++ b/test/numeric/test_regression_random_spikes.py @@ -35,7 +35,9 @@ def test_regression_random(): spike = spk.spike_distance_multi(spike_trains) spike_prof = spk.spike_profile_multi(spike_trains).avrg() - # spike_sync = spk.spike_sync_multi(spike_trains) + + spike_sync = spk.spike_sync_multi(spike_trains) + spike_sync_prof = spk.spike_sync_profile_multi(spike_trains).avrg() assert_almost_equal(isi, results_cSPIKY[i][0], decimal=14, err_msg="Index: %d, ISI" % i) @@ -47,6 +49,9 @@ def test_regression_random(): assert_almost_equal(spike_prof, results_cSPIKY[i][1], decimal=14, err_msg="Index: %d, SPIKE" % i) + assert_almost_equal(spike_sync, spike_sync_prof, decimal=14, + err_msg="Index: %d, SPIKE-Sync" % i) + def check_regression_dataset(spike_file="benchmark.mat", spikes_name="spikes", @@ -109,19 +114,25 @@ def check_single_spike_train_set(index): spike_train_data = spike_train_sets[index] spike_trains = [] + N = 0 for spikes in spike_train_data[0]: - print("Spikes:", spikes.flatten()) - spike_trains.append(spk.SpikeTrain(spikes.flatten(), 100.0)) + N += len(spikes.flatten()) + print("Spikes:", len(spikes.flatten())) + spikes_array = spikes.flatten() + if len(spikes_array > 0) and (spikes_array[-1] > 100.0): + spikes_array[-1] = 100.0 + spike_trains.append(spk.SpikeTrain(spikes_array, 100.0)) + print(spike_trains[-1].spikes) - print(spk.spike_distance_multi(spike_trains)) + print(N) - print(results_cSPIKY[index][1]) + print(spk.spike_sync_multi(spike_trains)) - print(spike_trains[1].spikes) + print(spk.spike_sync_profile_multi(spike_trains).integral()) if __name__ == "__main__": - test_regression_random() + # test_regression_random() # check_regression_dataset() - # check_single_spike_train_set(7633) + check_single_spike_train_set(4) diff --git a/test/test_distance.py b/test/test_distance.py index 083d8a3..fe09f34 100644 --- a/test/test_distance.py +++ b/test/test_distance.py @@ -217,6 +217,28 @@ def test_spike_sync(): assert_almost_equal(spk.spike_sync(spikes1, spikes2), 0.4, decimal=16) + spikes1 = SpikeTrain([1.0, 2.0, 4.0], 4.0) + spikes2 = SpikeTrain([3.8], 4.0) + spikes3 = SpikeTrain([3.9, ], 4.0) + + expected_x = np.array([0.0, 1.0, 2.0, 3.8, 4.0, 4.0]) + expected_y = np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0]) + + f = spk.spike_sync_profile(spikes1, spikes2) + + assert_array_almost_equal(f.x, expected_x, decimal=16) + assert_array_almost_equal(f.y, expected_y, decimal=16) + + f2 = spk.spike_sync_profile(spikes2, spikes3) + + i1 = f.integral() + i2 = f2.integral() + f.add(f2) + i12 = f.integral() + + assert_equal(i1[0]+i2[0], i12[0]) + assert_equal(i1[1]+i2[1], i12[1]) + def check_multi_profile(profile_func, profile_func_multi, dist_func_multi): # generate spike trains: @@ -313,6 +335,17 @@ def test_multi_spike_sync(): assert_equal(np.sum(f.y[1:-1]), 39932) assert_equal(np.sum(f.mp[1:-1]), 85554) + # example with 2 empty spike trains + sts = [] + sts.append(SpikeTrain([1, 9], [0, 10])) + sts.append(SpikeTrain([1, 3], [0, 10])) + sts.append(SpikeTrain([], [0, 10])) + sts.append(SpikeTrain([], [0, 10])) + + assert_almost_equal(spk.spike_sync_multi(sts), 1.0/6.0, decimal=15) + assert_almost_equal(spk.spike_sync_profile_multi(sts).avrg(), 1.0/6.0, + decimal=15) + def check_dist_matrix(dist_func, dist_matrix_func): # generate spike trains: diff --git a/test/test_save_load.py b/test/test_save_load.py new file mode 100644 index 0000000..33249b3 --- /dev/null +++ b/test/test_save_load.py @@ -0,0 +1,42 @@ +""" test_save_load.py + +Tests saving and loading of spike trains + +Copyright 2016, Mario Mulansky <mario.mulansky@gmx.net> + +Distributed under the BSD License + +""" + +from __future__ import print_function +from numpy.testing import assert_array_equal + +import tempfile +import os.path + +import pyspike as spk + + +def test_save_load(): + file_name = os.path.join(tempfile.mkdtemp(prefix='pyspike_'), + "save_load.txt") + + N = 10 + # generate some spike trains + spike_trains = [] + for n in range(N): + spike_trains.append(spk.generate_poisson_spikes(1.0, [0, 100])) + + # save them into txt file + spk.save_spike_trains_to_txt(spike_trains, file_name, precision=17) + + # load again + spike_trains_loaded = spk.load_spike_trains_from_txt(file_name, [0, 100]) + + for n in range(N): + assert_array_equal(spike_trains[n].spikes, + spike_trains_loaded[n].spikes) + + +if __name__ == "__main__": + test_save_load() |