further adjustments in spike sync

6 files changed, 49 insertions, 87 deletions

diff --git a/pyspike/SpikeTrain.py b/pyspike/SpikeTrain.py
index 041f897..89520c9 100644
--- a/pyspike/SpikeTrain.py
+++ b/pyspike/SpikeTrain.py
@@ -9,7 +9,7 @@ import numpy as np
 import collections
 
 
-class SpikeTrain:
+class SpikeTrain(object):
     """ Class representing spike trains for the PySpike Module."""
 
     def __init__(self, spike_times, interval):
diff --git a/pyspike/__init__.py b/pyspike/__init__.py
index 76e58a1..a5f9f0a 100644
--- a/pyspike/__init__.py
+++ b/pyspike/__init__.py
@@ -21,5 +21,5 @@ 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 add_auxiliary_spikes, load_spike_trains_from_txt, \
-    spike_train_from_string, merge_spike_trains, generate_poisson_spikes
+from spikes import load_spike_trains_from_txt, spike_train_from_string, \
+    merge_spike_trains, generate_poisson_spikes
diff --git a/pyspike/cython/cython_distance.pyx b/pyspike/cython/cython_distance.pyx
index 6d998b9..2841da8 100644
--- a/pyspike/cython/cython_distance.pyx
+++ b/pyspike/cython/cython_distance.pyx
@@ -337,10 +337,10 @@ def spike_distance_cython(double[:] t1, double[:] t2,
 # coincidence_python
 ############################################################
 cdef inline double get_tau(double[:] spikes1, double[:] spikes2,
-                           int i, int j, max_tau):
+                           int i, int j, double max_tau):
     cdef double m = 1E100   # some huge number
-    cdef int N1 = len(spikes1)-1
-    cdef int N2 = len(spikes2)-1
+    cdef int N1 = spikes1.shape[0]-1  # len(spikes1)-1
+    cdef int N2 = spikes2.shape[0]-1  # len(spikes2)-1
     if i < N1 and i > -1:
         m = fmin(m, spikes1[i+1]-spikes1[i])
     if j < N2 and j > -1:
diff --git a/pyspike/spike_sync.py b/pyspike/spike_sync.py
index bca6f73..8ddd32c 100644
--- a/pyspike/spike_sync.py
+++ b/pyspike/spike_sync.py
@@ -109,10 +109,10 @@ def spike_sync_profile_multi(spike_trains, indices=None, max_tau=None):
 
     """
     prof_func = partial(spike_sync_profile, max_tau=max_tau)
-    average_dist, M = _generic_profile_multi(spike_trains, prof_func,
+    average_prof, M = _generic_profile_multi(spike_trains, prof_func,
                                              indices)
     # average_dist.mul_scalar(1.0/M)  # no normalization here!
-    return average_dist
+    return average_prof
 
 
 ############################################################
@@ -122,7 +122,7 @@ def spike_sync_multi(spike_trains, indices=None, interval=None, max_tau=None):
     """ Computes the multi-variate spike synchronization value for a set of
     spike trains.
 
-    :param spike_trains: list of spike trains
+    :param spike_trains: list of :class:`pyspike.SpikeTrain`
     :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
diff --git a/pyspike/spikes.py b/pyspike/spikes.py
index 9d7d6f4..128873d 100644
--- a/pyspike/spikes.py
+++ b/pyspike/spikes.py
@@ -8,82 +8,46 @@ Distributed under the BSD License
 """
 
 import numpy as np
-
-
-############################################################
-# add_auxiliary_spikes
-############################################################
-def add_auxiliary_spikes(spike_train, time_interval):
-    """ Adds spikes at the beginning and end of the given time interval.
-
-    :param spike_train: ordered array of spike times
-    :param time_interval: A pair (T_start, T_end) of values representing the
-                          start and end time of the spike train measurement or
-                          a single value representing the end time, the T_start
-                          is then assuemd as 0. Auxiliary spikes will be added
-                          to the spike train at the beginning and end of this
-                          interval, if they are not yet present.
-    :type time_interval: pair of doubles or double
-    :returns: spike train with additional spikes at T_start and T_end.
-
-    """
-    try:
-        T_start = time_interval[0]
-        T_end = time_interval[1]
-    except:
-        T_start = 0
-        T_end = time_interval
-
-    assert spike_train[0] >= T_start, \
-        "Spike train has events before the given start time"
-    assert spike_train[-1] <= T_end, \
-        "Spike train has events after the given end time"
-    if spike_train[0] != T_start:
-        spike_train = np.insert(spike_train, 0, T_start)
-    if spike_train[-1] != T_end:
-        spike_train = np.append(spike_train, T_end)
-    return spike_train
+from pyspike import SpikeTrain
 
 
 ############################################################
 # spike_train_from_string
 ############################################################
-def spike_train_from_string(s, sep=' ', is_sorted=False):
-    """ Converts a string of times into an array of spike times.
+def spike_train_from_string(s, interval, sep=' ', is_sorted=False):
+    """ Converts a string of times into a  :class:`pyspike.SpikeTrain`.
 
-    :param s: the string with (ordered) spike times
+    :param s: the string with (ordered) spike times.
+    :param interval: interval defining the edges of the spike train.
+    Given as a pair of floats (T0, T1) or a single float T1, where T0=0 is
+    assumed.
     :param sep: The separator between the time numbers, default=' '.
     :param is_sorted: if True, the spike times are not sorted after loading,
                       if False, spike times are sorted with `np.sort`
-    :returns: array of spike times
+    :returns:  :class:`pyspike.SpikeTrain`
     """
     if not(is_sorted):
-        return np.sort(np.fromstring(s, sep=sep))
+        return SpikeTrain(np.sort(np.fromstring(s, sep=sep)), interval)
     else:
-        return np.fromstring(s, sep=sep)
+        return SpikeTrain(np.fromstring(s, sep=sep), interval)
 
 
 ############################################################
 # load_spike_trains_txt
 ############################################################
-def load_spike_trains_from_txt(file_name, time_interval=None,
+def load_spike_trains_from_txt(file_name, interval=None,
                                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
-    neglected. The `time_interval` represents the start and the end of the
-    spike trains and it is used to add auxiliary spikes at the beginning and
-    end of each spike train. However, if `time_interval == None`, no auxiliary
-    spikes are added, but note that the Spike and ISI distance both require
-    auxiliary spikes.
+    neglected. The `interval` represents the start and the end of the
+    spike trains.
 
     :param file_name: The name of the text file.
-    :param time_interval: A pair (T_start, T_end) of values representing the
-                          start and end time of the spike train measurement
-                          or a single value representing the end time, the
-                          T_start is then assuemd as 0. Auxiliary spikes will
-                          be added to the spike train at the beginning and end
-                          of this interval.
+    :param interval: A pair (T_start, T_end) of values representing the
+                     start and end time of the spike train measurement
+                     or a single value representing the end time, the
+                     T_start is then assuemd as 0.
     :param separator: The character used to seprate the values in the text file
     :param comment: Lines starting with this character are ignored.
     :param sort: If true, the spike times are order via `np.sort`, default=True
@@ -94,9 +58,8 @@ 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, is_sorted)
-            if time_interval is not None:  # add auxil. spikes if times given
-                spike_train = add_auxiliary_spikes(spike_train, time_interval)
+            spike_train = spike_train_from_string(line, interval,
+                                                  separator, is_sorted)
             spike_trains.append(spike_train)
     return spike_trains
 
@@ -111,14 +74,14 @@ def merge_spike_trains(spike_trains):
     :returns: spike train with the merged spike times
     """
     # get the lengths of the spike trains
-    lens = np.array([len(st) for st in spike_trains])
+    lens = np.array([len(st.spikes) for st in spike_trains])
     merged_spikes = np.empty(np.sum(lens))
     index = 0                             # the index for merged_spikes
     indices = np.zeros_like(lens)         # indices of the spike trains
     index_list = np.arange(len(indices))  # indices of indices of spike trains
                                           # that have not yet reached the end
     # list of the possible events in the spike trains
-    vals = [spike_trains[i][indices[i]] for i in index_list]
+    vals = [spike_trains[i].spikes[indices[i]] for i in index_list]
     while len(index_list) > 0:
         i = np.argmin(vals)             # the next spike is the minimum
         merged_spikes[index] = vals[i]  # put it to the merged spike train
@@ -127,33 +90,34 @@ def merge_spike_trains(spike_trains):
         indices[i] += 1                 # next index for the chosen spike train
         if indices[i] >= lens[i]:       # remove spike train index if ended
             index_list = index_list[index_list != i]
-        vals = [spike_trains[n][indices[n]] for n in index_list]
-    return merged_spikes
+        vals = [spike_trains[n].spikes[indices[n]] for n in index_list]
+    return SpikeTrain(merged_spikes, [spike_trains[0].t_start,
+                                      spike_trains[0].t_end])
 
 
 ############################################################
 # generate_poisson_spikes
 ############################################################
-def generate_poisson_spikes(rate, time_interval, add_aux_spikes=True):
+def generate_poisson_spikes(rate, interval):
     """ Generates a Poisson spike train with the given rate in the given time
     interval
 
     :param rate: The rate of the spike trains
-    :param time_interval: A pair (T_start, T_end) of values representing the
-                          start and end time of the spike train measurement or
-                          a single value representing the end time, the T_start
-                          is then assuemd as 0. Auxiliary spikes will be added
-                          to the spike train at the beginning and end of this
-                          interval, if they are not yet present.
-    :type time_interval: pair of doubles or double
-    :returns: Poisson spike train
+    :param interval: A pair (T_start, T_end) of values representing the
+                     start and end time of the spike train measurement or
+                     a single value representing the end time, the T_start
+                     is then assuemd as 0. Auxiliary spikes will be added
+                     to the spike train at the beginning and end of this
+                     interval, if they are not yet present.
+    :type interval: pair of doubles or double
+    :returns: Poisson spike train as a :class:`pyspike.SpikeTrain`
     """
     try:
-        T_start = time_interval[0]
-        T_end = time_interval[1]
+        T_start = interval[0]
+        T_end = interval[1]
     except:
         T_start = 0
-        T_end = time_interval
+        T_end = interval
     # roughly how many spikes are required to fill the interval
     N = max(1, int(1.2 * rate * (T_end-T_start)))
     N_append = max(1, int(0.1 * rate * (T_end-T_start)))
@@ -165,7 +129,4 @@ def generate_poisson_spikes(rate, time_interval, add_aux_spikes=True):
                               np.random.exponential(1.0/rate, N_append))
     spikes = T_start + np.cumsum(intervals)
     spikes = spikes[spikes < T_end]
-    if add_aux_spikes:
-        return add_auxiliary_spikes(spikes, time_interval)
-    else:
-        return spikes
+    return SpikeTrain(spikes, interval)
diff --git a/test/test_distance.py b/test/test_distance.py
index dbb72f1..0fff840 100644
--- a/test/test_distance.py
+++ b/test/test_distance.py
@@ -250,7 +250,8 @@ def test_multi_spike_sync():
 
     # multivariate regression test
     spike_trains = spk.load_spike_trains_from_txt("test/SPIKE_Sync_Test.txt",
-                                                  time_interval=(0, 4000))
+                                                  interval=(0, 4000))
+    print(spike_trains[0].spikes)
     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)
@@ -339,4 +340,4 @@ if __name__ == "__main__":
     test_spike_sync()
     test_multi_isi()
     test_multi_spike()
-    test_multi_spike_sync()
+    # test_multi_spike_sync()