Repository: NeuralEnsemble/elephant

Issue:
[Feautre]  build with py 3.8
Hi team ---

Does not work with py 3.8

Cheers,

Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index acceba1..987bbec 100644
--- a/.gitignore
+++ b/.gitignore
@@ -66,7 +66,7 @@ cover
 ######################
 .directory
 .gdb_history
-.DS_Store?
+.DS_Store
 ehthumbs.db
 Icon?
 Thumbs.db
@@ -75,3 +75,5 @@ Thumbs.db
 ###################################
 # ignored folder for fast prototyping
 ignored/
+
+.ipynb_checkpoints/
diff --git a/.travis.yml b/.travis.yml
index e383982..864f58c 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -21,14 +21,6 @@ matrix:
       python: 3.5
       env: DISTRIB="pip"
 
-    - name: "conda 3.6"
-      python: 3.6
-      env: DISTRIB="conda"
-
-    - name: "pip 3.6"
-      python: 3.6
-      env: DISTRIB="pip"
-
     - name: "pip 3.6 requirements-extras"
       python: 3.6
       env: DISTRIB="pip"
@@ -41,6 +33,14 @@ matrix:
       python: 3.7
       env: DISTRIB="conda"
 
+    - name: "conda 3.8"
+      python: 3.8
+      env: DISTRIB="conda"
+
+    - name: "pip 3.8"
+      python: 3.8
+      env: DISTRIB="pip"
+
   exclude:
     - name: "pip 3.6 requirements-extras min version"
       # excluded due to unrelated unittest assertWarn bug error
diff --git a/elephant/conversion.py b/elephant/conversion.py
index ff9ac35..8c2b981 100644
--- a/elephant/conversion.py
+++ b/elephant/conversion.py
@@ -832,7 +832,6 @@ class BinnedSpikeTrain(object):
     def remove_stored_array(self):
         """
         Unlinks the matrix with counted time points from memory.
-
         """
         self._mat_u = None
 
diff --git a/elephant/spade.py b/elephant/spade.py
index 0f83fd4..2d4e3e7 100644
--- a/elephant/spade.py
+++ b/elephant/spade.py
@@ -56,31 +56,34 @@ plt.show()
 :copyright: Copyright 2017 by the Elephant team, see `doc/authors.rst`.
 :license: BSD, see LICENSE.txt for details.
 """
-import numpy as np
-from scipy import sparse
+import time
+import warnings
 import operator
-import neo
-import elephant.spike_train_surrogates as surr
-import elephant.conversion as conv
 from itertools import chain, combinations
 from functools import reduce
-import time
+from collections import defaultdict
+
+import numpy as np
+import neo
 import quantities as pq
-import warnings
+from scipy import sparse
+
+import elephant.spike_train_surrogates as surr
+import elephant.conversion as conv
 from elephant.spade_src import fast_fca
-from collections import defaultdict
 
 warnings.simplefilter('once', UserWarning)
 
-
 try:
     from mpi4py import MPI  # for parallelized routines
+
     HAVE_MPI = True
 except ImportError:  # pragma: no cover
     HAVE_MPI = False
 
 try:
     from elephant.spade_src import fim
+
     HAVE_FIM = True
 except ImportError:  # pragma: no cover
     HAVE_FIM = False
@@ -90,7 +93,7 @@ def spade(data, binsize, winlen, min_spikes=2, min_occ=2, max_spikes=None,
           max_occ=None, min_neu=1, n_subsets=0, delta=0, epsilon=0,
           stability_thresh=None, n_surr=0, dither=15 * pq.ms, spectrum='#',
           alpha=1, stat_corr='fdr', psr_param=None, output_format='concepts'):
-    """
+    r"""
     Perform the SPADE [1,2] analysis for the parallel spike trains given in the
     input. The data are discretized with a temporal resolution equal binsize
     in a sliding window of winlen*binsize milliseconds.
@@ -190,11 +193,11 @@ def spade(data, binsize, winlen, min_spikes=2, min_occ=2, max_spikes=None,
     psr_param: None or list of int
         This list contains parameters used in the pattern spectrum filtering:
             psr_param[0]: correction parameter for subset filtering
-                (see parameter h of pattern_set_reduction()).
+                (see h_subset_filtering in pattern_set_reduction()).
             psr_param[1]: correction parameter for superset filtering
-                (see parameter k of pattern_set_reduction()).
+                (see k_superset_filtering in pattern_set_reduction()).
             psr_param[2]: correction parameter for covered-spikes criterion
-                (see parameter l for pattern_set_reduction()).
+                (see l_covered_spikes in pattern_set_reduction()).
     output_format: str
         distinguish the format of the output (see Returns). Can assume values
         'concepts' and 'patterns'.
@@ -285,41 +288,30 @@ def spade(data, binsize, winlen, min_spikes=2, min_occ=2, max_spikes=None,
     Frontiers in Computational Neuroscience, 11.
     """
     if HAVE_MPI:  # pragma: no cover
-        comm = MPI.COMM_WORLD   # create MPI communicator
+        comm = MPI.COMM_WORLD  # create MPI communicator
         rank = comm.Get_rank()  # get rank of current MPI task
     else:
         rank = 0
-    output = {}
-    # Decide whether compute pvalue spectrum
-    if n_surr > 0:
-        # Compute pvalue spectrum
-        time_pvalue_spectrum = time.time()
-        pv_spec = pvalue_spectrum(data, binsize, winlen, dither=dither,
-                                  n_surr=n_surr, min_spikes=min_spikes,
-                                  min_occ=min_occ, max_spikes=max_spikes,
-                                  max_occ=max_occ, min_neu=min_neu,
-                                  spectrum=spectrum)
-        time_pvalue_spectrum = time.time() - time_pvalue_spectrum
-        print("Time for pvalue spectrum computation: {}".format(
-            time_pvalue_spectrum))
-        # Storing pvalue spectrum
-        output['pvalue_spectrum'] = pv_spec
-    elif 0 < alpha < 1:
-        warnings.warn('0<alpha<1 but p-value spectrum has not been '
-                      'computed (n_surr==0)')
+
+    if output_format not in ['concepts', 'patterns']:
+        raise ValueError("The output_format value has to be"
+                         "'patterns' or 'concepts'")
+
     time_mining = time.time()
-    # Decide if compute the approximated stability
-    if n_subsets > 0:
+    if rank == 0 or n_subsets > 0:
         # Mine the data for extraction of concepts
         concepts, rel_matrix = concepts_mining(data, binsize, winlen,
                                                min_spikes=min_spikes,
                                                min_occ=min_occ,
-                                               min_neu=min_neu,
                                                max_spikes=max_spikes,
                                                max_occ=max_occ,
+                                               min_neu=min_neu,
                                                report='a')
         time_mining = time.time() - time_mining
         print("Time for data mining: {}".format(time_mining))
+
+    # Decide if compute the approximated stability
+    if n_subsets > 0:
         # Computing the approximated stability of all the concepts
         time_stability = time.time()
         concepts = approximate_stability(concepts, rel_matrix, n_subsets,
@@ -333,69 +325,72 @@ def spade(data, binsize, winlen, min_spikes=2, min_occ=2, max_spikes=None,
     elif stability_thresh is not None:
         warnings.warn('Stability_thresh not None but stability has not been '
                       'computed (n_subsets==0)')
-    elif rank == 0:
-        # Mine the data for extraction of concepts
-        concepts, rel_matrix = concepts_mining(data, binsize, winlen,
-                                               min_spikes=min_spikes,
-                                               min_occ=min_occ,
-                                               max_spikes=max_spikes,
-                                               max_occ=max_occ,
-                                               min_neu=min_neu,
-                                               report='a')
-        time_mining = time.time() - time_mining
-        print("Time for data mining: {}".format(time_mining))
-    if rank == 0:
-        # Decide whether filter concepts with psf
-        if n_surr > 0:
-            if len(pv_spec) == 0:
-                ns_sgnt = []
-            else:
-                # Computing non-significant entries of the spectrum applying
-                # the statistical correction
-                ns_sgnt = test_signature_significance(pv_spec, alpha,
-                                                      corr=stat_corr,
-                                                      report='non_significant',
-                                                      spectrum=spectrum)
-            # Storing non-significant entries of the pvalue spectrum
-            output['non_sgnf_sgnt'] = ns_sgnt
-            # Filter concepts with pvalue spectrum (psf)
-            if len(ns_sgnt) > 0:
-                concepts = list(filter(
-                    lambda c: _pattern_spectrum_filter(
-                        c, ns_sgnt, spectrum, winlen), concepts))
-                # Decide whether to filter concepts using psr
-                if psr_param is not None:
-                    # Filter using conditional tests (psr)
-                    concepts = pattern_set_reduction(concepts, ns_sgnt,
-                                                     winlen=winlen,
-                                                     h=psr_param[0],
-                                                     k=psr_param[1],
-                                                     l=psr_param[2],
-                                                     min_spikes=min_spikes,
-                                                     min_occ=min_occ)
-
-        # Storing patterns
-        if output_format == 'patterns':
-            # If the p-value spectra was not computed, is set to an empty list
-            if n_surr == 0:
-                pv_spec = None
-            # Transforming concepts to dictionary containing pattern's infos
-            output['patterns'] = concept_output_to_patterns(concepts,
-                                                            winlen, binsize,
-                                                            pv_spec, spectrum,
-                                                            data[0].t_start)
-        elif output_format == 'concepts':
-            output['patterns'] = concepts
-        else:
-            raise ValueError(
-                "The output_format value has to be one between"
-                "'patterns' and 'concepts'")
-        return output
+
+    output = {}
+    pv_spec = None  # initialize pv_spec to None
+    # Decide whether compute pvalue spectrum
+    if n_surr > 0:
+        # Compute pvalue spectrum
+        time_pvalue_spectrum = time.time()
+        pv_spec = pvalue_spectrum(data, binsize, winlen, dither=dither,
+                                  n_surr=n_surr, min_spikes=min_spikes,
+                                  min_occ=min_occ, max_spikes=max_spikes,
+                                  max_occ=max_occ, min_neu=min_neu,
+                                  spectrum=spectrum)
+        time_pvalue_spectrum = time.time() - time_pvalue_spectrum
+        print("Time for pvalue spectrum computation: {}".format(
+            time_pvalue_spectrum))
+        # Storing pvalue spectrum
+        output['pvalue_spectrum'] = pv_spec
+    elif 0 < alpha < 1:
+        warnings.warn('0<alpha<1 but p-value spectrum has not been '
+                      'computed (n_surr==0)')
+
     # rank!=0 returning None
-    else:
+    if rank != 0:
         warnings.warn('Returning None because executed on a process != 0')
         return None
 
+    # Initialize non-significant signatures as empty list:
+    ns_sgnt = []
+    # Decide whether filter concepts with psf
+    if n_surr > 0:
+        if len(pv_spec) > 0:
+            # Computing non-significant entries of the spectrum applying
+            # the statistical correction
+            ns_sgnt = test_signature_significance(pv_spec, alpha,
+                                                  corr=stat_corr,
+                                                  report='non_significant',
+                                                  spectrum=spectrum)
+        # Storing non-significant entries of the pvalue spectrum
+        output['non_sgnf_sgnt'] = ns_sgnt
+    # Filter concepts with pvalue spectrum (psf)
+    if len(ns_sgnt) > 0:
+        concepts = list(filter(
+            lambda c: _pattern_spectrum_filter(
+                c, ns_sgnt, spectrum, winlen), concepts))
+        # Decide whether to filter concepts using psr
+        if psr_param is not None:
+            # Filter using conditional tests (psr)
+            concepts = pattern_set_reduction(concepts, ns_sgnt,
+                                             winlen=winlen,
+                                             h_subset_filtering=psr_param[0],
+                                             k_superset_filtering=psr_param[1],
+                                             l_covered_spikes=psr_param[2],
+                                             min_spikes=min_spikes,
+                                             min_occ=min_occ)  # nopep8
+    # Storing patterns for ouput format concepts
+    if output_format == 'concepts':
+        output['patterns'] = concepts
+        return output
+
+    # Transforming concepts to dictionary containing pattern's infos
+    output['patterns'] = concept_output_to_patterns(concepts,
+                                                    winlen, binsize,
+                                                    pv_spec, spectrum,
+                                                    data[0].t_start)
+    return output
+
 
 def concepts_mining(data, binsize, winlen, min_spikes=2, min_occ=2,
                     max_spikes=None, max_occ=None, min_neu=1, report='a'):
@@ -599,7 +594,6 @@ def _build_context(binary_matrix, winlen, only_windows_with_first_spike=True):
     # Shape of the rel_matrix:
     # (num of window positions,
     # num of bins in one window * number of neurons)
-    num_windows = window_indices.shape[0]
     rel_matrix = sparse.coo_matrix(
         (np.ones((len(windows_col)), dtype=bool),
          (windows_row, windows_col)),
@@ -714,73 +708,72 @@ def _fpgrowth(transactions, min_c=2, min_z=2, max_z=None,
             and first spike of the pattern, number of patterns)
 
     """
+    if min_neu < 1:
+        raise AttributeError('min_neu must be an integer >=1')
     # By default, set the maximum pattern size to the number of spiketrains
     if max_z is None:
         max_z = np.max((np.max([len(tr) for tr in transactions]), min_z + 1))
     # By default set maximum number of data to number of bins
     if max_c is None:
         max_c = len(transactions)
-    if min_neu >= 1:
-        # Initializing outputs
-        concepts = []
-        if report == '#':
-            spec_matrix = np.zeros((max_z + 1, max_c + 1))
-        if report == '3d#':
-            spec_matrix = np.zeros((max_z + 1, max_c + 1, winlen))
-        spectrum = []
-        # Mining the data with fpgrowth algorithm
-        if np.unique(transactions, return_counts=True)[1][0] == len(
-                transactions):
-            fpgrowth_output = [(tuple(transactions[0]), len(transactions))]
-        else:
-            fpgrowth_output = fim.fpgrowth(
-                tracts=transactions,
-                target=target,
-                supp=-min_c,
-                zmin=min_z,
-                zmax=max_z,
-                report='a',
-                algo='s')
-        # Applying min/max conditions and computing extent (window positions)
-        fpgrowth_output = list(filter(
-            lambda c: _fpgrowth_filter(
-                c, winlen, max_c, min_neu), fpgrowth_output))
-        # filter out subsets of patterns that are found as a side-effect
-        # of using the moving window strategy
-        fpgrowth_output = _filter_for_moving_window_subsets(
-            fpgrowth_output, winlen)
-        for (intent, supp) in fpgrowth_output:
-            if report == 'a':
-                if rel_matrix is not None:
-                    # Computing the extent of the concept (patterns
-                    # occurrences), checking in rel_matrix in which windows
-                    # the intent occurred
-                    extent = tuple(np.where(
-                        np.all(rel_matrix[:, intent], axis=1) == 1)[0])
-                concepts.append((intent, extent))
-            # Computing 2d spectrum
-            elif report == '#':
-                spec_matrix[len(intent) - 1, supp - 1] += 1
-            # Computing 3d spectrum
-            elif report == '3d#':
-                spec_matrix[len(intent) - 1, supp - 1, max(
-                    np.array(intent) % winlen)] += 1
-        del fpgrowth_output
-        if report == 'a':
-            return concepts
-        else:
-            if report == '#':
-                for (z, c) in np.transpose(np.where(spec_matrix != 0)):
-                    spectrum.append((z + 1, c + 1, int(spec_matrix[z, c])))
-
-            elif report == '3d#':
-                for (z, c, l) in np.transpose(np.where(spec_matrix != 0)):
-                    spectrum.append(
-                        (z + 1, c + 1, l, int(spec_matrix[z, c, l])))
-            del spec_matrix
-            return spectrum
+
+    # Initializing outputs
+    concepts = []
+    if report == '#':
+        spec_matrix = np.zeros((max_z + 1, max_c + 1))
+    if report == '3d#':
+        spec_matrix = np.zeros((max_z + 1, max_c + 1, winlen))
+    spectrum = []
+    # Mining the data with fpgrowth algorithm
+    if np.unique(transactions, return_counts=True)[1][0] == len(
+            transactions):
+        fpgrowth_output = [(tuple(transactions[0]), len(transactions))]
     else:
-        raise AttributeError('min_neu must be an integer >=1')
+        fpgrowth_output = fim.fpgrowth(
+            tracts=transactions,
+            target=target,
+            supp=-min_c,
+            zmin=min_z,
+            zmax=max_z,
+            report='a',
+            algo='s')
+    # Applying min/max conditions and computing extent (window positions)
+    fpgrowth_output = list(filter(
+        lambda c: _fpgrowth_filter(
+            c, winlen, max_c, min_neu), fpgrowth_output))
+    # filter out subsets of patterns that are found as a side-effect
+    # of using the moving window strategy
+    fpgrowth_output = _filter_for_moving_window_subsets(
+        fpgrowth_output, winlen)
+    for (intent, supp) in fpgrowth_output:
+        if report == 'a':
+            if rel_matrix is not None:
+                # Computing the extent of the concept (patterns
+                # occurrences), checking in rel_matrix in which windows
+                # the intent occurred
+                extent = tuple(np.where(
+                    np.all(rel_matrix[:, intent], axis=1) == 1)[0])
+            concepts.append((intent, extent))
+        # Computing 2d spectrum
+        elif report == '#':
+            spec_matrix[len(intent) - 1, supp - 1] += 1
+        # Computing 3d spectrum
+        elif report == '3d#':
+            spec_matrix[len(intent) - 1, supp - 1, max(
+                np.array(intent) % winlen)] += 1
+    del fpgrowth_output
+    if report == 'a':
+        return concepts
+
+    if report == '#':
+        for (z, c) in np.transpose(np.where(spec_matrix != 0)):
+            spectrum.append((z + 1, c + 1, int(spec_matrix[z, c])))
+    elif report == '3d#':
+        for (z, c, l) in np.transpose(np.where(spec_matrix != 0)):
+            spectrum.append(
+                (z + 1, c + 1, l, int(spec_matrix[z, c, l])))
+    del spec_matrix
+    return spectrum
 
 
 def _fpgrowth_filter(concept, winlen, max_c, min_neu):
@@ -792,7 +785,7 @@ def _fpgrowth_filter(concept, winlen, max_c, min_neu):
     keep_concepts = len(
         np.unique(np.array(
             concept[0]) // winlen)) >= min_neu and concept[1] <= max_c and min(
-                np.array(concept[0]) % winlen) == 0
+        np.array(concept[0]) % winlen) == 0
     return keep_concepts
 
 
@@ -977,19 +970,19 @@ def _fast_fca(context, min_c=2, min_z=2, max_z=None,
                 np.array(intent) % winlen)] += 1
     if report == 'a':
         return concepts
-    else:
-        del concepts
-        # returning spectrum
-        if report == '#':
-            for (z, c) in np.transpose(np.where(spec_matrix != 0)):
-                spectrum.append((z + 1, c + 1, int(spec_matrix[z, c])))
 
-        if report == '3d#':
-            for (z, c, l) in np.transpose(np.where(spec_matrix != 0)):
-                spectrum.append(
-                    (z + 1, c + 1, l, int(spec_matrix[z, c, l])))
-        del spec_matrix
-        return spectrum
+    del concepts
+    # returning spectrum
+    if report == '#':
+        for (z, c) in np.transpose(np.where(spec_matrix != 0)):
+            spectrum.append((z + 1, c + 1, int(spec_matrix[z, c])))
+
+    if report == '3d#':
+        for (z, c, l) in np.transpose(np.where(spec_matrix != 0)):
+            spectrum.append(
+                (z + 1, c + 1, l, int(spec_matrix[z, c, l])))
+    del spec_matrix
+    return spectrum
 
 
 def _fca_filter(concept, winlen, min_c, min_z, max_c, max_z, min_neu):
@@ -1001,8 +994,8 @@ def _fca_filter(concept, winlen, min_c, min_z, max_c, max_z, min_neu):
     extent = tuple(concept.extent)
     keep_concepts = len(intent) >= min_z and len(extent) >= min_c and len(
         intent) <= max_z and len(extent) <= max_c and len(
-            np.unique(np.array(intent) // winlen)) >= min_neu and min(
-                np.array(intent) % winlen) == 0
+        np.unique(np.array(intent) // winlen)) >= min_neu and min(
+        np.array(intent) % winlen) == 0
     return keep_concepts
 
 
@@ -1086,7 +1079,7 @@ def pvalue_spectrum(data, binsize, winlen, dither, n_surr, min_spikes=2,
     """
     # Initializing variables for parallel computing
     if HAVE_MPI:  # pragma: no cover
-        comm = MPI.COMM_WORLD   # create MPI communicator
+        comm = MPI.COMM_WORLD  # create MPI communicator
         rank = comm.Get_rank()  # get rank of current MPI task
         size = comm.Get_size()  # get tot number of MPI tasks
     else:
@@ -1190,7 +1183,7 @@ def _fdr(pvalues, alpha):
     for i, pv in enumerate(pvs_sorted):  # For each PV, from the largest on
         k = m - i
         if pv <= alpha * (k * 1. / m):  # continue if PV > fdr-threshold
-            break                          # otherwise stop
+            break  # otherwise stop
     # this applies, when loop is not stopped due to significant pvalue
     else:
         k = 0
@@ -1278,6 +1271,15 @@ def test_signature_significance(pvalue_spectrum, alpha, corr='',
     # If alpha == 1 all signatures are significant
     if alpha == 1:
         return []
+
+    if spectrum not in ['#', '3d#']:
+        raise AttributeError("spectrum must be either '#' or '3d#', "
+                             "got {} instead".format(spectrum))
+    if report not in ['spectrum', 'significant', 'non_significant']:
+        raise AttributeError("report must be either 'spectrum'," +
+                             "  'significant' or 'non_significant'," +
+                             "got {} instead".format(report))
+
     x_array = np.array(pvalue_spectrum)
     # Compute significance...
     if corr in ['', 'no']:  # ...without statistical correction
@@ -1291,39 +1293,31 @@ def test_signature_significance(pvalue_spectrum, alpha, corr='',
     else:
         raise AttributeError(
             "Parameter corr must be either ''('no'), 'b'('bonf'), 'f'('fdr')" +
-            " or 'hb'('holm_bonf'), but not '"+str(corr)+"'.")
+            " or 'hb'('holm_bonf'), but not '" + str(corr) + "'.")
     # Return the specified results:
     if spectrum == '#':
         if report == 'spectrum':
             return [(size, supp, test)
                     for (size, supp, pv), test in zip(pvalue_spectrum, tests)]
-        elif report == 'significant':
+        if report == 'significant':
             return [(size, supp) for ((size, supp, pv), test)
                     in zip(pvalue_spectrum, tests) if test]
-        elif report == 'non_significant':
-            return [
-                (size, supp) for ((size, supp, pv), test) in zip(
-                    pvalue_spectrum, tests) if not test]
-        raise AttributeError("report must be either 'spectrum'," +
-                             "  'significant' or 'non_significant'," +
-                             "got {} instead".format(report))
-    elif spectrum == '3d#':
-        if report == 'spectrum':
-            return [(size, supp, l, test)
-                    for (size, supp, l, pv), test in zip(
-                    pvalue_spectrum, tests)]
-        elif report == 'significant':
-            return [(size, supp, l) for ((size, supp, l, pv), test)
-                    in zip(pvalue_spectrum, tests) if test]
-        elif report == 'non_significant':
-            return [
-                (size, supp, l) for ((size, supp, l, pv), test) in zip(
-                    pvalue_spectrum, tests) if not test]
-        raise AttributeError("report must be either 'spectrum'," +
-                             "  'significant' or 'non_significant'," +
-                             "got {} instead".format(report))
-    raise AttributeError("spectrum must be either '#' or '3d#', "
-                         "got {} instead".format(spectrum))
+        # report == 'non_significant'
+        return [(size, supp)
+                for ((size, supp, pv), test) in zip(pvalue_spectrum, tests)
+                if not test]
+
+    # spectrum == '3d#'
+    if report == 'spectrum':
+        return [(size, supp, l, test)
+                for (size, supp, l, pv), test in zip(pvalue_spectrum, tests)]
+    if report == 'significant':
+        return [(size, supp, l) for ((size, supp, l, pv), test)
+                in zip(pvalue_spectrum, tests) if test]
+    # report == 'non_significant'
+    return [(size, supp, l)
+            for ((size, supp, l, pv), test) in zip(pvalue_spectrum, tests)
+            if not test]
 
 
 def _pattern_spectrum_filter(concept, ns_signature, spectrum, winlen):
@@ -1345,7 +1339,7 @@ def _pattern_spectrum_filter(concept, ns_signature, spectrum, winlen):
 
 
 def approximate_stability(concepts, rel_matrix, n_subsets, delta=0, epsilon=0):
-    """
+    r"""
     Approximate the stability of concepts. Uses the algorithm described
     in Babin, Kuznetsov (2012): Approximating Concept Stability
 
@@ -1409,7 +1403,7 @@ def approximate_stability(concepts, rel_matrix, n_subsets, delta=0, epsilon=0):
 
     """
     if HAVE_MPI:  # pragma: no cover
-        comm = MPI.COMM_WORLD   # create MPI communicator
+        comm = MPI.COMM_WORLD  # create MPI communicator
         rank = comm.Get_rank()  # get rank of current MPI task
         size = comm.Get_size()  # get tot number of MPI tasks
     else:
@@ -1419,7 +1413,7 @@ def approximate_stability(concepts, rel_matrix, n_subsets, delta=0, epsilon=0):
         raise AttributeError('n_subsets has to be >=0 or delta + epsilon > 0')
     if len(concepts) == 0:
         return []
-    elif len(concepts) <= size:
+    if len(concepts) <= size:
         rank_idx = [0] * (size + 1) + [len(concepts)]
     else:
         rank_idx = list(
@@ -1429,154 +1423,86 @@ def approximate_stability(concepts, rel_matrix, n_subsets, delta=0, epsilon=0):
     # Calculate optimal n
     if delta + epsilon > 0 and n_subsets == 0:
         n_subsets = np.log(2. / delta) / (2 * epsilon ** 2) + 1
-    output = []
+
     if rank == 0:
-        for concept in concepts[
-                rank_idx[rank]:rank_idx[rank + 1]] + concepts[
-                rank_idx[-2]:rank_idx[-1]]:
-            stab_ext = 0.0
-            stab_int = 0.0
-            intent = np.array(list(concept[0]))
-            extent = np.array(list(concept[1]))
-            r_unique_ext = set()
-            r_unique_int = set()
-            excluded_subset = []
-            # Calculate all subsets if n is larger than the power set of
-            # the extent
-            if n_subsets > 2 ** len(extent):
-                subsets_ext = chain.from_iterable(
-                    combinations(extent, r) for r in range(
-                        len(extent) + 1))
-                for s in subsets_ext:
-                    if any(
-                            [set(s).issubset(se) for se in excluded_subset]):
-                        continue
-                    if _closure_probability_extensional(
-                            intent, s, rel_matrix):
-                        stab_ext += 1
-                    else:
-                        excluded_subset.append(s)
-            else:
-                for _ in range(n_subsets):
-                    subset_ext = extent[
-                        _give_random_idx(r_unique_ext, len(extent))]
-                    if any([
-                        set(subset_ext).issubset(se) for
-                            se in excluded_subset]):
-                        continue
-                    if _closure_probability_extensional(
-                            intent, subset_ext, rel_matrix):
-                        stab_ext += 1
-                    else:
-                        excluded_subset.append(subset_ext)
-            stab_ext /= min(n_subsets, 2 ** len(extent))
-            excluded_subset = []
-            # Calculate all subsets if n is larger than the power set of
-            # the extent
-            if n_subsets > 2 ** len(intent):
-                subsets_int = chain.from_iterable(
-                    combinations(intent, r) for r in range(
-                        len(intent) + 1))
-                for s in subsets_int:
-                    if any(
-                            [set(s).issubset(se) for se in excluded_subset]):
-                        continue
-                    if _closure_probability_intensional(
-                            extent, s, rel_matrix):
-                        stab_int += 1
-                    else:
-                        excluded_subset.append(s)
-            else:
-                for _ in range(n_subsets):
-                    subset_int = intent[
-                        _give_random_idx(r_unique_int, len(intent))]
-                    if any([
-                        set(subset_int).issubset(se) for
-                            se in excluded_subset]):
-                        continue
-                    if _closure_probability_intensional(
-                            extent, subset_int, rel_matrix):
-                        stab_int += 1
-                    else:
-                        excluded_subset.append(subset_int)
-            stab_int /= min(n_subsets, 2 ** len(intent))
-            output.append((intent, extent, stab_int, stab_ext))
-    else:  # pragma: no cover
-        for concept in concepts[rank_idx[rank]:rank_idx[rank + 1]]:
-            stab_ext = 0.0
-            stab_int = 0.0
-            intent = np.array(list(concept[0]))
-            extent = np.array(list(concept[1]))
-            r_unique_ext = set()
-            r_unique_int = set()
-            excluded_subset = []
-            # Calculate all subsets if n is larger than the power set of
-            # the extent
-            if n_subsets > 2 ** len(extent):
-                subsets_ext = chain.from_iterable(
-                    combinations(extent, r) for r in range(
-                        len(extent) + 1))
-                for s in subsets_ext:
-                    if any(
-                            [set(s).issubset(se) for se in excluded_subset]):
-                        continue
-                    if _closure_probability_extensional(
-                            intent, s, rel_matrix):
-                        stab_ext += 1
-                    else:
-                        excluded_subset.append(s)
-            else:
-                for _ in range(n_subsets):
-                    subset_ext = extent[
-                        _give_random_idx(r_unique_ext, len(extent))]
-                    if any([
-                        set(subset_ext).issubset(se) for
-                            se in excluded_subset]):
-                        continue
-                    if _closure_probability_extensional(
-                            intent, subset_ext, rel_matrix):
-                        stab_ext += 1
-                    else:
-                        excluded_subset.append(subset_ext)
-            stab_ext /= min(n_subsets, 2 ** len(extent))
-            excluded_subset = []
-            # Calculate all subsets if n is larger than the power set of
-            # the extent
-            if n_subsets > 2 ** len(intent):
-                subsets_int = chain.from_iterable(
-                    combinations(intent, r) for r in range(
-                        len(intent) + 1))
-                for s in subsets_int:
-                    if any(
-                            [set(s).issubset(se) for se in excluded_subset]):
-                        continue
-                    if _closure_probability_intensional(
-                            extent, s, rel_matrix):
-                        stab_int += 1
-                    else:
-                        excluded_subset.append(s)
-            else:
-                for _ in range(n_subsets):
-                    subset_int = intent[
-                        _give_random_idx(r_unique_int, len(intent))]
-                    if any([
-                        set(subset_int).issubset(se) for
-                            se in excluded_subset]):
-                        continue
-                    if _closure_probability_intensional(
-                            extent, subset_int, rel_matrix):
-                        stab_int += 1
-                    else:
-                        excluded_subset.append(subset_int)
-            stab_int /= min(n_subsets, 2 ** len(intent))
-            output.append((intent, extent, stab_int, stab_ext))
+        concepts_on_partition = concepts[rank_idx[rank]:rank_idx[rank + 1]] + \
+            concepts[rank_idx[-2]:rank_idx[-1]]
+    else:
+        concepts_on_partition = concepts[rank_idx[rank]:rank_idx[rank + 1]]
 
-    if rank != 0:  # pragma: no cover
-        comm.send(output, dest=0)
-    if rank == 0:  # pragma: no cover
-        for i in range(1, size):
-            recv_list = comm.recv(source=i)
-            output.extend(recv_list)
+    output = []
+    for concept in concepts_on_partition:
+        stab_ext = 0.0
+        stab_int = 0.0
+        intent = np.array(list(concept[0]))
+        extent = np.array(list(concept[1]))
+        r_unique_ext = set()
+        r_unique_int = set()
+        excluded_subset = []
+        # Calculate all subsets if n is larger than the power set of
+        # the extent
+        if n_subsets > 2 ** len(extent):
+            subsets_ext = chain.from_iterable(
+                combinations(extent, r) for r in range(
+                    len(extent) + 1))
+            for s in subsets_ext:
+                if any(
+                        [set(s).issubset(se) for se in excluded_subset]):
+                    continue
+                if _closure_probability_extensional(
+                        intent, s, rel_matrix):
+                    stab_ext += 1
+                else:
+                    excluded_subset.append(s)
+        else:
+            for _ in range(n_subsets):
+                subset_ext = extent[
+                    _give_random_idx(r_unique_ext, len(extent))]
+                if any([set(subset_ext).issubset(se)
+                        for se in excluded_subset]):
+                    continue
+                if _closure_probability_extensional(
+                        intent, subset_ext, rel_matrix):
+                    stab_ext += 1
+                else:
+                    excluded_subset.append(subset_ext)
+        stab_ext /= min(n_subsets, 2 ** len(extent))
+        excluded_subset = []
+        # Calculate all subsets if n is larger than the power set of
+        # the extent
+        if n_subsets > 2 ** len(intent):
+            subsets_int = chain.from_iterable(
+                combinations(intent, r) for r in range(
+                    len(intent) + 1))
+            for s in subsets_int:
+                if any(
+                        [set(s).issubset(se) for se in excluded_subset]):
+                    continue
+                if _closure_probability_intensional(
+                        extent, s, rel_matrix):
+                    stab_int += 1
+                else:
+                    excluded_subset.append(s)
+        else:
+            for _ in range(n_subsets):
+                subset_int = intent[
+                    _give_random_idx(r_unique_int, len(intent))]
+                if any([set(subset_int).issubset(se) for
+                        se in excluded_subset]):
+                    continue
+                if _closure_probability_intensional(
+                        extent, subset_int, rel_matrix):
+                    stab_int += 1
+                else:
+                    excluded_subset.append(subset_int)
+        stab_int /= min(n_subsets, 2 ** len(intent))
+        output.append((intent, extent, stab_int, stab_ext))
+
+    if size != 1:
+        recv_list = comm.gather(output, root=0)
+        if rank == 0:
+            for i in range(1, len(recv_list)):
+                output.extend(recv_list[i])
 
     return output
 
@@ -1637,8 +1563,7 @@ def _give_random_idx(r_unique, n):
     """ asd """
 
     r = np.random.randint(n,
-                          size=np.random.randint(low=1,
-                                                 high=n))
+                          size=np.random.randint(low=1, high=n))
     r_tuple = tuple(r)
     if r_tuple not in r_unique:
         r_unique.add(r_tuple)
@@ -1646,9 +1571,10 @@ def _give_random_idx(r_unique, n):
     return _give_random_idx(r_unique, n)
 
 
-def pattern_set_reduction(concepts, excluded, winlen, h=0, k=0, l=0,
+def pattern_set_reduction(concepts, excluded, winlen, h_subset_filtering=0,
+                          k_superset_filtering=0, l_covered_spikes=0,
                           min_spikes=2, min_occ=2):
-    """
+    r"""
     Takes a list concepts and performs pattern set reduction (PSR).
 
     PSR determines which patterns in concepts_psf are statistically significant
@@ -1682,13 +1608,13 @@ def pattern_set_reduction(concepts, excluded, winlen, h=0, k=0, l=0,
         List of concepts, each consisting in its intent and extent
     excluded: list
         A list of non-significant pattern signatures (z, c) (see above).
-    h: int
+    h_subset_filtering: int
         Correction parameter for subset filtering (see above).
         Defaults: 0
-    k: int
+    k_superset_filtering: int
         Correction parameter for superset filtering (see above).
         Default: 0
-    l: int ]
+    l_covered_spikes: int
         Correction parameter for covered-spikes criterion (see above).
         Default: 0
     min_size: int
@@ -1743,21 +1669,21 @@ def pattern_set_reduction(concepts, excluded, winlen, h=0, k=0, l=0,
                     continue
                 if not selected[id2]:
                     continue
-                if set(conc1_new).issuperset(conc2) and count2\
-                        - count1 + h < min_occ:
+                if set(conc1_new).issuperset(conc2) and count2 \
+                        - count1 + h_subset_filtering < min_occ:
                     selected[id2] = False
                     break
-                if set(conc2).issuperset(conc1_new) and count1\
-                        - count2 + h < min_occ:
+                if set(conc2).issuperset(conc1_new) and count1 \
+                        - count2 + h_subset_filtering < min_occ:
                     selected[id1] = False
                     break
                 if len(excluded) == 0:
                     break
                 # Test the case conc1 is a superset of conc2
                 if set(conc1_new).issuperset(conc2):
-                    supp_diff = count2 - count1 + h
+                    supp_diff = count2 - count1 + h_subset_filtering
                     size1, size2 = len(conc1_new), len(conc2)
-                    size_diff = size1 - size2 + k
+                    size_diff = size1 - size2 + k_superset_filtering
                     # 2d spectrum case
                     if len(excluded[0]) == 2:
                         # Determine whether the subset (conc2)
@@ -1788,23 +1714,22 @@ def pattern_set_reduction(concepts, excluded, winlen, h=0, k=0, l=0,
                     if reject_sub and not reject_sup:
                         selected[id2] = False
                         break
-                    elif reject_sup and not reject_sub:
+                    if reject_sup and not reject_sub:
                         selected[id1] = False
                         break
-                    elif reject_sub and reject_sup:
-                        if (size1 - l) * count1 >= (size2 - l) * count2:
+                    if reject_sub and reject_sup:
+                        if (size1 - l_covered_spikes) * \
+                                count1 >= (size2 - l_covered_spikes) * count2:
                             selected[id2] = False
                             break
-                        else:
-                            selected[id1] = False
-                            break
+                        selected[id1] = False
+                        break
                     # if both sets are significant given the other, keep both
-                    else:
-                        continue
+                    continue
                 elif set(conc2).issuperset(conc1_new):
-                    supp_diff = count1 - count2 + h
+                    supp_diff = count1 - count2 + h_subset_filtering
                     size1, size2 = len(conc1_new), len(conc2)
-                    size_diff = size2 - size1 + k
+                    size_diff = size2 - size1 + k_superset_filtering
                     # 2d spectrum case
                     if len(excluded[0]) == 2:
                         # Determine whether the subset (conc2) should be
@@ -1816,7 +1741,7 @@ def pattern_set_reduction(concepts, excluded, winlen, h=0, k=0, l=0,
                         reject_sup = (size_diff, count1) in excluded \
                             or size_diff < min_spikes
                     # 3d spectrum case
-                    if len(excluded[0]) == 3:
+                    elif len(excluded[0]) == 3:
                         # Determine whether the subset (conc2) should be
                         # rejected according to the test for excess occurrences
                         len_sub = max(
@@ -1834,63 +1759,66 @@ def pattern_set_reduction(concepts, excluded, winlen, h=0, k=0, l=0,
                     if reject_sub and not reject_sup:
                         selected[id1] = False
                         break
-                    elif reject_sup and not reject_sub:
+                    if reject_sup and not reject_sub:
                         selected[id2] = False
                         break
-                    elif reject_sub and reject_sup:
-                        if (size1 - l) * count1 >= (size2 - l) * count2:
+                    if reject_sub and reject_sup:
+                        if (size1 - l_covered_spikes) * \
+                                count1 >= (size2 - l_covered_spikes) * count2:
                             selected[id2] = False
                             break
-                        else:
-                            selected[id1] = False
-                            break
+                        selected[id1] = False
+                        break
                     # if both sets are significant given the other, keep both
-                    else:
-                        continue
+                    continue
                 else:
                     size1, size2 = len(conc1_new), len(conc2)
                     inter_size = len(set(conc1_new) & set(conc2))
                     # 2d spectrum case
                     if len(excluded[0]) == 2:
-                        reject_1 = (
-                            size1 - inter_size + k,
-                            count1) in \
-                            excluded or size1 - inter_size + k < min_spikes
-                        reject_2 = (
-                            size2 - inter_size + k, count2) in excluded or \
-                            size2 - inter_size + k < min_spikes
+                        reject_1 = (size1 - inter_size + k_superset_filtering,
+                                    count1) in excluded
+                        reject_1 = reject_1 or \
+                            (size1 - inter_size + k_superset_filtering)\
+                            < min_spikes
+                        reject_2 = (size2 - inter_size + k_superset_filtering,
+                                    count2) in excluded
+                        reject_2 = reject_2 or \
+                            (size1 - inter_size + k_superset_filtering) \
+                            < min_spikes
                     # 3d spectrum case
                     if len(excluded[0]) == 3:
                         len_1 = max(
                             np.abs(np.diff(np.array(conc1_new) % winlen)))
                         len_2 = max(
                             np.abs(np.diff(np.array(conc2) % winlen)))
-                        reject_1 = (
-                            size1 - inter_size + k, count1,
-                            len_1) in excluded or \
-                            size1 - inter_size + k < min_spikes
-                        reject_2 = (
-                            size2 - inter_size + k, count2, len_2) \
-                            in excluded or \
-                            size2 - inter_size + k < min_spikes
+                        reject_1 = (size1 - inter_size + k_superset_filtering,
+                                    count1, len_1) in excluded
+                        reject_1 = reject_1 or \
+                            (size1 - inter_size + k_superset_filtering) \
+                            < min_spikes
+                        reject_2 = (size2 - inter_size + k_superset_filtering,
+                                    count2, len_2) in excluded
+                        reject_2 = reject_2 or \
+                            (size1 - inter_size + k_superset_filtering) \
+                            < min_spikes
 
                     # Reject accordingly:
                     if reject_2 and not reject_1:
                         selected[id2] = False
                         break
-                    elif reject_1 and not reject_2:
+                    if reject_1 and not reject_2:
                         selected[id1] = False
                         break
-                    elif reject_1 and reject_2:
-                        if (size1 - l) * count1 >= (size2 - l) * count2:
+                    if reject_1 and reject_2:
+                        if (size1 - l_covered_spikes) * \
+                                count1 >= (size2 - l_covered_spikes) * count2:
                             selected[id2] = False
                             break
-                        else:
-                            selected[id1] = False
-                            break
+                        selected[id1] = False
+                        break
                     # if both sets are significant given the other, keep both
-                    else:
-                        continue
+                    continue
 
     # Return the selected concepts
     return [p for i, p in enumerate(concepts) if selected[i]]
diff --git a/elephant/spike_train_generation.py b/elephant/spike_train_generation.py
index 9820698..82dc13b 100644
--- a/elephant/spike_train_generation.py
+++ b/elephant/spike_train_generation.py
@@ -1159,7 +1159,7 @@ def homogeneous_poisson_process_with_refr_period(rate,
 
     # Check that the number of spikes drawn from the Poisson distribution,
     # can fit in the duration regarding the refractory period.
-    spike_count = min(spike_count, math.ceil(duration / refr_period))
+    spike_count = min(spike_count, math.ceil(duration / refr_period_mag))
 
     # Due to the refractory period the effective space in where the spikes
     # can be placed is shortened by (spike_count - 1) times the
diff --git a/elephant/spike_train_surrogates.py b/elephant/spike_train_surrogates.py
index 6eb445c..a8e19dc 100644
--- a/elephant/spike_train_surrogates.py
+++ b/elephant/spike_train_surrogates.py
@@ -26,6 +26,10 @@ the original data:
     keeping the first spike time fixed and generating the others from the
     new sequence of ISIs. Keeps spike count and ISIs, flattens the firing rate
     profile
+* joint_isi_dithering:
+    calculate the Joint-ISI distribution and moves spike according to the
+    probability distribution, that results from a fixed sum of ISI_before
+    and the ISI_afterwards. For further details see [1].
 
 [1] Louis et al (2010) Surrogate Spike Train Generation Through Dithering in
     Operational Time. Front Comput Neurosci. 2010; 4: 127.
@@ -38,8 +42,12 @@ Original implementation by: Emiliano Torre [e.torre@fz-juelich.de]
 import numpy as np
 import quantities as pq
 import neo
+from scipy.ndimage import gaussian_filter
+import random
+
 try:
     import elephant.statistics as es
+
     isi = es.isi
 except ImportError:
     from .statistics import isi  # Convenience when in elephant working dir.
@@ -100,7 +108,8 @@ def dither_spikes(spiketrain, dither, n=1, decimals=None, edges=True):
         818.84446913]) * ms, [0.0 ms, 1000.0 ms])>,
      <SpikeTrain(array([ 111.36693058,  235.15750163,  618.87388515,
         786.1807108 ]) * ms, [0.0 ms, 1000.0 ms])>]
-    >>> print dither_spikes(st, dither = 20*pq.ms, decimals=0)   # doctest: +SKIP
+    >>> print dither_spikes(st, dither = 20*pq.ms,
+                            decimals=0)   # doctest: +SKIP
     [<SpikeTrain(array([  81.,  242.,  595.,  799.]) * ms,
         [0.0 ms, 1000.0 ms])>]
     """
@@ -119,12 +128,12 @@ def dither_spikes(spiketrain, dither, n=1, decimals=None, edges=True):
         # Move all spikes outside [spiketrain.t_start, spiketrain.t_stop] to
         # the range's ends
         surr = np.minimum(np.maximum(surr.simplified.magnitude,
-            spiketrain.t_start.simplified.magnitude),
-            spiketrain.t_stop.simplified.magnitude) * pq.s
+                                     spiketrain.t_start.simplified.magnitude),
+                          spiketrain.t_stop.simplified.magnitude) * pq.s
     else:
         # Leave out all spikes outside [spiketrain.t_start, spiketrain.t_stop]
         tstart, tstop = spiketrain.t_start.simplified.magnitude, \
-                        spiketrain.t_stop.simplified.magnitude
+            spiketrain.t_stop.simplified.magnitude
         surr = [np.sort(s[np.all([s >= tstart, s < tstop], axis=0)]) * pq.s
                 for s in surr.simplified.magnitude]
 
@@ -194,7 +203,8 @@ def randomise_spikes(spiketrain, n=1, decimals=None):
         sts = sts.round(decimals)
 
     # Convert the Quantity array to a list of SpikeTrains, and return them
-    return [neo.SpikeTrain(np.sort(st), t_start=spiketrain.t_start, t_stop=spiketrain.t_stop)
+    return [neo.SpikeTrain(np.sort(st), t_start=spiketrain.t_start,
+                           t_stop=spiketrain.t_stop)
             for st in sts]
 
 
@@ -255,7 +265,7 @@ def shuffle_isis(spiketrain, n=1, decimals=None):
         # Create list of surrogate spike trains by random ISI permutation
         sts = []
         for i in range(n):
-            surr_times = np.cumsum(np.random.permutation(ISIs)) *\
+            surr_times = np.cumsum(np.random.permutation(ISIs)) * \
                 spiketrain.units + spiketrain.t_start
             sts.append(neo.SpikeTrain(
                 surr_times, t_start=spiketrain.t_start,
@@ -327,7 +337,8 @@ def dither_spike_train(spiketrain, shift, n=1, decimals=None, edges=True):
         792.89084054]) * ms, [0.0 ms, 1000.0 ms])>,
      <SpikeTrain(array([  84.61079043,  234.61079043,  584.61079043,
         784.61079043]) * ms, [0.0 ms, 1000.0 ms])>]
-    >>> print dither_spike_train(st, shift = 20*pq.ms, decimals=0)   # doctest: +SKIP
+    >>> print dither_spike_train(st, shift = 20*pq.ms,
+                                 decimals=0)   # doctest: +SKIP
     [<SpikeTrain(array([  82.,  232.,  582.,  782.]) * ms,
         [0.0 ms, 1000.0 ms])>]
     """
@@ -336,8 +347,8 @@ def dither_spike_train(spiketrain, shift, n=1, decimals=None, edges=True):
     data = spiketrain.view(pq.Quantity)
 
     # Main: generate the surrogates by spike train shifting
-    surr = data.reshape((1, len(data))) + 2 * shift * \
-        np.random.random_sample((n, 1)) - shift
+    surr = data.reshape(
+        (1, len(data))) + 2 * shift * np.random.random_sample((n, 1)) - shift
 
     # Round the surrogate data to decimal position, if requested
     if decimals is not None:
@@ -347,12 +358,12 @@ def dither_spike_train(spiketrain, shift, n=1, decimals=None, edges=True):
         # Move all spikes outside [spiketrain.t_start, spiketrain.t_stop] to
         # the range's ends
         surr = np.minimum(np.maximum(surr.simplified.magnitude,
-            spiketrain.t_start.simplified.magnitude),
-            spiketrain.t_stop.simplified.magnitude) * pq.s
+                                     spiketrain.t_start.simplified.magnitude),
+                          spiketrain.t_stop.simplified.magnitude) * pq.s
     else:
         # Leave out all spikes outside [spiketrain.t_start, spiketrain.t_stop]
         tstart, tstop = spiketrain.t_start.simplified.magnitude, \
-                        spiketrain.t_stop.simplified.magnitude
+            spiketrain.t_stop.simplified.magnitude
         surr = [np.sort(s[np.all([s >= tstart, s < tstop], axis=0)]) * pq.s
                 for s in surr.simplified.magnitude]
 
@@ -372,7 +383,7 @@ def jitter_spikes(spiketrain, binsize, n=1):
 
     The surrogates retain the :attr:`t_start and :attr:`t_stop` of the
     :attr:`spike train`. Note that within each time bin the surrogate
-    `neo.SpikeTrain` objects are locally poissonian (the inter-spike-interval
+    `neo.SpikeTrain` objects are locally Poissonian (the inter-spike-interval
     are exponentially distributed).
 
     Parameters
@@ -441,7 +452,7 @@ def jitter_spikes(spiketrain, binsize, n=1):
     offsets = start_dl + np.array([bin_edges[bin_id] for bin_id in bin_ids])
     dilats = np.array([bin_sizes_dl[bin_id] for bin_id in bin_ids])
 
-    # Compute each surrogate by dilatating and shifting each spike s in the
+    # Compute each surrogate by dilating and shifting each spike s in the
     # poisson 0-1 spike trains to dilat * s + offset. Attach time unit again
     surr = np.sort(surr_poiss01 * dilats + offsets, axis=1) * std_unit
 
@@ -450,9 +461,451 @@ def jitter_spikes(spiketrain, binsize, n=1):
             for s in surr]
 
 
+class JointISI(object):
+    """
+    The class :class:`JointISI` is implemented for Joint-ISI dithering
+    as a continuation of the ideas of Louis et al. (2010) and Gerstein (2004).
+
+    When creating an class instance all necessary preprocessing steps are done,
+    to use the method dithering().
+
+    Attributes
+    ----------
+    spiketrain: neo.SpikeTrain
+        For this spiketrain the surrogates will be created
+
+    dither: pq.Quantity
+        This quantity describes the maximum displacement of a spike, when
+        method is 'window'. It is also used for the uniform dithering for
+        the spikes, which are outside the regime in the Joint-ISI
+        histogram, where Joint-ISI dithering is applicable.
+        Default: 15.*pq.ms
+    truncation_limit: pq.Quantity
+        The Joint-ISI distribution is as such defined on a range for ISI_i and
+        ISI_(i+1) from 0 to inf. Since this is computationally not feasible,
+        the Joint-ISI distribution is truncated for high ISI. The Joint-ISI
+        histogram is calculated for ISI_i, ISI_(i+1) from 0 to
+        truncation_limit.
+        Default: 100*pq.ms
+    num_bins: int
+        The size of the joint-ISI-distribution will be num_bins*num_bins/2.
+        Default: 100
+    sigma: pq.Quantity
+        The standard deviation of the Gaussian kernel, with which
+        the data is convoluted.
+        Default: 2.*pq.ms
+    alternate: boolean
+        If alternate == True: then first all even and then all odd spikes are
+        dithered. Else: in ascending order from the first to the last spike,
+        all spikes are moved.
+        Default: True.
+    use_sqrt: boolean
+        if use_sqrt == True a sqrt is applied to the joint-ISI histogram,
+        following Gerstein et al. 2004
+        Default: False
+    method: string
+        if 'window': the spike movement is limited to the parameter dither.
+        if 'fast': the spike can move in all the range between the previous
+            spike and the subsequent spike. This is computationally much faster
+            and thus is called 'fast'.
+        Default: 'fast'
+    cutoff: boolean
+        if True then the Filtering of the Joint-ISI histogram is
+        limited to the lower side by the minimal ISI.
+        This can be necessary, if in the data there is a certain refractory
+        period, which would be destroyed by the convolution with the
+        2d-Gaussian function.
+        Default: True
+    refr_period:
+        Since this dither-method should conserve the refractory period,
+        It is internally calculated as the minimum of the value given here
+        and the smallest ISI in the spiketrain.
+        Default: 4.*pq.ms
+    Methods
+    ----------
+    dithering()
+        Returns a list of dithered spiketrains.
+    """
+
+    # The min number of spikes, required for dithering.
+    # Otherwise, the original spiketrain is copied N times.
+    MIN_SPIKES = 3
+
+    def __init__(self,
+                 spiketrain,
+                 dither=15. * pq.ms,
+                 truncation_limit=100. * pq.ms,
+                 num_bins=100,
+                 sigma=2. * pq.ms,
+                 alternate=True,
+                 use_sqrt=False,
+                 method='fast',
+                 cutoff=True,
+                 refr_period=4. * pq.ms
+                 ):
+        """
+        Parameters
+        ----------
+        spiketrain: neo.SpikeTrain
+            Input spiketrain to create surrogates from.
+        dither: pq.Quantity, optional
+            This quantity describes the maximum displacement of a spike, when
+            method is 'window'. It is also used for the uniform dithering for
+            the spikes, which are outside the regime in the Joint-ISI
+            histogram, where Joint-ISI dithering is applicable.
+            Default: 15.*pq.ms
+        truncation_limit: pq.Quantity, optional
+            The Joint-ISI distribution of :math:`(ISI_i, ISI_{i+1})` is defined
+            within the range `[0, inf]`. Since this is computationally not
+            feasible, the Joint-ISI distribution is truncated for high ISI.
+            The Joint-ISI histogram is calculated for
+            :math:`(ISI_i, ISI_{i+1})` from 0 to `truncation_limit`.
+            Default: 100*pq.ms
+        num_bins: int, optional
+            The size of the joint-ISI-distribution will be
+            `num_bins*num_bins/2`.
+            Default: 100
+        sigma: pq.Quantity, optional
+            The standard deviation of the Gaussian kernel, with which
+            the data is convolved.
+            Default: 2.*pq.ms
+        alternate: boolean, optional
+            If `alternate` is True, then all even spikes are dithered followed
+            by all odd spikes. Otherwise, the spikes are dithered in ascending
+            order.
+            Default: True.
+        use_sqrt: boolean, optional
+            If `use_sqrt` is True, the joint-ISI histogram is preprocessed with
+            a square root (following Gerstein et al. 2004).
+            Default: False
+        method: string, optional
+            * 'window': the spike movement is limited to the parameter `dither`
+            * 'fast': the spike can move in the whole range between the
+                previous and subsequent spikes (computationally efficient).
+            Default: 'fast'
+        cutoff: boolean, optional
+            If set to True, then the Filtering of the Joint-ISI histogram is
+            limited to the lower side by the minimal ISI.
+            This can be necessary, if in the data there is a certain refractory
+            period, which will be destroyed by the convolution with the
+            2d-Gaussian function.
+            Default: True
+        refr_period: pq.Quantity, optional
+            Defines the refractory period of the dithered `spiketrain` unless
+            the smallest ISI of the `spiketrain` is lower than this value.
+            Default: 4.*pq.ms
+        """
+        self.spiketrain = spiketrain
+        self.truncation_limit = self.get_magnitude(truncation_limit)
+        self.num_bins = num_bins
+
+        self.dither = self.get_magnitude(dither)
+
+        self.sigma = self.get_magnitude(sigma)
+        self.alternate = alternate
+
+        if method not in ['fast', 'window']:
+            raise ValueError("The method can either be 'fast' or 'window', "
+                             "but not '{}'".format(method))
+        self.method = method
+
+        refr_period = self.get_magnitude(refr_period)
+        if not self.too_less_spikes:
+            minimal_isi = np.min(self.isi)
+            refr_period = min(refr_period, minimal_isi)
+        self.refr_period = refr_period
+
+        self.cutoff = cutoff
+        self.use_sqrt = use_sqrt
+        self._jisih_cumulatives = None
+
+        self.max_change_index = self.isi_to_index(self.dither)
+        self.max_change_isi = self.index_to_isi(self.max_change_index)
+
+    def get_magnitude(self, x):
+        """
+        Parameters
+        ----------
+        x: pq.Quantity or float
+
+        Returns
+        -------
+        float
+            The magnitude of `x`, rescaled to the units of the input
+            `spiketrain`.
+        """
+        if isinstance(x, pq.Quantity):
+            return x.rescale(self.unit).magnitude
+        return x
+
+    @property
+    def too_less_spikes(self):
+        """
+        This is a check if the spiketrain has enough spikes to evaluate the
+        joint-ISI histogram.
+
+        Returns
+        -------
+        bool
+        """
+        return len(self.spiketrain) < self.MIN_SPIKES
+
+    @property
+    def unit(self):
+        return self.spiketrain.units
+
+    @property
+    def isi(self):
+        if self.too_less_spikes:
+            return None
+        return isi(self.spiketrain.magnitude)
+
+    @property
+    def bin_width(self):
+        return self.truncation_limit / self.num_bins
+
+    def isi_to_index(self, inter_spike_interval):
+        """
+        A function that gives for each ISI the corresponding index in the
+        Joint-ISI distribution.
+
+        Parameters
+        ----------
+        inter_spike_interval: np.ndarray or float
+            An array of ISIs or a single ISI.
+
+        Returns
+        -------
+        np.ndarray or int:
+            The corresponding index for each ISI.
+        """
+        return np.floor(inter_spike_interval / self.bin_width).astype(int)
+
+    def index_to_isi(self, isi_index):
+        """
+        Maps `isi_index` back to the original ISI.
+
+        Parameters
+        ----------
+        isi_index: np.ndarray or int
+            The index of ISI.
+
+        Returns
+        -------
+        np.ndarray or float:
+            The corresponding ISI for each index.
+        """
+        return (isi_index + 0.5) * self.bin_width
+
+    def joint_isi_histogram(self):
+        if self.too_less_spikes:
+            return None
+        isis = self.isi
+        joint_isi_histogram = np.histogram2d(
+            isis[:-1], isis[1:],
+            bins=[self.num_bins, self.num_bins],
+            range=[[0., self.truncation_limit],
+                   [0., self.truncation_limit]])[0]
+
+        if self.use_sqrt:
+            joint_isi_histogram = np.sqrt(joint_isi_histogram)
+
+        if self.sigma:
+            if self.cutoff:
+                start_index = self.isi_to_index(self.refr_period)
+                joint_isi_histogram[
+                    start_index:, start_index:] = gaussian_filter(
+                    joint_isi_histogram[start_index:, start_index:],
+                    sigma=self.sigma / self.bin_width)
+                joint_isi_histogram[:start_index, :] = 0
+                joint_isi_histogram[:, :start_index] = 0
+            else:
+                joint_isi_histogram = gaussian_filter(
+                    joint_isi_histogram, sigma=self.sigma / self.bin_width)
+        return joint_isi_histogram
+
+    @staticmethod
+    def normalize_cumulative_distribution(array):
+        """
+        This function normalizes parts of a cumulative distribution function
+        to be a cumulative distribution function again.
+
+        Parameters
+        ----------
+        array: np.ndarray
+            A monotonously increasing array as a part of an unnormalized
+             cumulative distribution function.
+        Returns
+        -------
+        np.ndarray
+            Monotonously increasing array from 0 to 1.
+        """
+        if array[-1] - array[0] > 0.:
+            return (array - array[0]) / (array[-1] - array[0])
+        return np.zeros_like(array)
+
+    def dithering(self, n_surr=1):
+        """
+        Implementation of Joint-ISI-dithering for spiketrains that pass the
+        threshold of the dense rate, if not a uniform dithered spiketrain is
+        given back. The implementation continued the ideas of Louis et al.
+        (2010) and Gerstein (2004).
+
+        Parameters
+        ----------
+        n_surr: int
+            The number of dithered spiketrains to be returned.
+            Default: 1
+
+        Returns
+        ----------
+        dithered_sts: list
+            List of spiketrains, that are dithered versions of the given
+            `spiketrain`
+        """
+        if self.too_less_spikes:
+            return [self.spiketrain] * n_surr
+
+        # Checks, whether the preprocessing is already done.
+        if self._jisih_cumulatives is None:
+            self._determine_cumulative_functions()
+
+        dithered_sts = []
+        isi_to_dither = self.isi
+        for surr_number in range(n_surr):
+            dithered_isi = self._get_dithered_isi(isi_to_dither)
+
+            dithered_st = self.spiketrain[0].magnitude + \
+                np.r_[0., np.cumsum(dithered_isi)]
+            dithered_st = neo.SpikeTrain(dithered_st * self.unit,
+                                         t_start=self.spiketrain.t_start,
+                                         t_stop=self.spiketrain.t_stop)
+            dithered_sts.append(dithered_st)
+        return dithered_sts
+
+    def _determine_cumulative_functions(self):
+        rotated_jisih = np.rot90(self.joint_isi_histogram())
+
+        if self.method == 'fast':
+            self._jisih_cumulatives = []
+            for double_index in range(self.num_bins):
+                # Taking anti-diagonals of the original joint-ISI histogram
+                diagonal = np.diagonal(
+                    rotated_jisih, offset=-self.num_bins + double_index + 1)
+                jisih_cum = self.normalize_cumulative_distribution(
+                    np.cumsum(diagonal))
+                self._jisih_cumulatives.append(jisih_cum)
+            self._jisih_cumulatives = np.array(self._jisih_cumulatives)
+        else:
+            self._jisih_cumulatives = self._window_cumulatives(rotated_jisih)
+
+    def _window_cumulatives(self, rotated_jisih):
+        jisih_diag_cums = self._window_diagonal_cumulatives(rotated_jisih)
+        jisih_cumulatives = np.zeros(
+            (self.num_bins, self.num_bins,
+             2 * self.max_change_index + 1))
+        for curr_isi_id in range(self.num_bins):
+            for next_isi_id in range(self.num_bins - curr_isi_id):
+                double_index = next_isi_id + curr_isi_id
+                cum_slice = jisih_diag_cums[
+                    double_index,
+                    curr_isi_id: curr_isi_id + 2 * self.max_change_index + 1]
+
+                normalized_cum = self.normalize_cumulative_distribution(
+                    cum_slice)
+                jisih_cumulatives[curr_isi_id][next_isi_id] = normalized_cum
+        return jisih_cumulatives
+
+    def _window_diagonal_cumulatives(self, rotated_jisih):
+        # An element of the first axis is defined as the sum of indices
+        # for previous and subsequent ISI.
+
+        jisih_diag_cums = np.zeros((self.num_bins,
+                                    self.num_bins
+                                    + 2 * self.max_change_index))
+
+        # double_index corresponds to the sum of the indices for the previous
+        # and the subsequent ISI.
+        for double_index in range(self.num_bins):
+            cum_diag = np.cumsum(np.diagonal(rotated_jisih,
+                                             - self.num_bins
+                                             + double_index + 1))
+
+            right_padding = jisih_diag_cums.shape[1] - \
+                len(cum_diag) - self.max_change_index
+
+            jisih_diag_cums[double_index] = np.pad(
+                cum_diag,
+                pad_width=(self.max_change_index, right_padding),
+                mode='constant',
+                constant_values=(cum_diag[0], cum_diag[-1])
+            )
+
+        return jisih_diag_cums
+
+    def _get_dithered_isi(self, isi_to_dither):
+        dithered_isi = np.copy(isi_to_dither)
+        # if alternate is true, a sampling_rhythm of 2 means that we have two
+        # partitions of spikes, the "odd" and the "even" spikes and first
+        # dither the "even" ones and then the "odd" ones.
+        # if alternate is false, we just go dither from the first to the last
+        # spike, which corresponds to a sampling_rhythm of 1.
+        sampling_rhythm = self.alternate + 1
+        number_of_isis = len(dithered_isi)
+
+        for start in range(sampling_rhythm):
+            dithered_isi_indices = self.isi_to_index(dithered_isi)
+            for i in range(start, number_of_isis - 1,
+                           sampling_rhythm):
+                step = self._get_dithering_step(
+                    dithered_isi,
+                    dithered_isi_indices,
+                    i)
+                dithered_isi[i] += step
+                dithered_isi[i + 1] -= step
+
+        return dithered_isi
+
+    def _get_dithering_step(self,
+                            dithered_isi,
+                            dithered_isi_indices,
+                            i):
+        curr_isi_id = dithered_isi_indices[i]
+        next_isi_id = dithered_isi_indices[i + 1]
+        double_index = curr_isi_id + next_isi_id
+        if double_index < self.num_bins:
+            if self.method == 'fast':
+                cum_dist_func = self._jisih_cumulatives[
+                    double_index]
+                compare_isi = self.index_to_isi(curr_isi_id)
+            else:
+                cum_dist_func = self._jisih_cumulatives[
+                    curr_isi_id][next_isi_id]
+                compare_isi = self.max_change_isi
+
+            if cum_dist_func[-1] > 0.:
+                # when the method is 'fast', new_isi_id is where the current
+                # ISI id should go to.
+                new_isi_id = np.searchsorted(cum_dist_func, random.random())
+                step = self.index_to_isi(new_isi_id) - compare_isi
+                return step
+
+        return self._uniform_dither_not_jisi_movable_spikes(
+            dithered_isi[i],
+            dithered_isi[i + 1])
+
+    def _uniform_dither_not_jisi_movable_spikes(self,
+                                                curr_isi,
+                                                next_isi):
+        left_dither = min(curr_isi - self.refr_period, self.dither)
+        right_dither = min(next_isi - self.refr_period, self.dither)
+        step = random.random() * (right_dither + left_dither) - left_dither
+        return step
+
+
 def surrogates(
-        spiketrain, n=1, surr_method='dither_spike_train', dt=None, decimals=None,
-        edges=True):
+        spiketrain, n=1, surr_method='dither_spike_train', dt=None,
+        decimals=None, edges=True):
     """
     Generates surrogates of a :attr:`spiketrain` by a desired generation
     method.
@@ -478,6 +931,7 @@ def surrogates(
         * 'jitter_spikes': see surrogates.jitter_spikes() [dt needed]
         * 'randomise_spikes': see surrogates.randomise_spikes()
         * 'shuffle_isis': see surrogates.shuffle_isis()
+        * 'joint_isi_dithering': see surrogates.joint_isi_dithering()
         Default: 'dither_spike_train'
     dt : quantities.Quantity, optional
         For methods shifting spike times randomly around their original time
@@ -509,7 +963,8 @@ def surrogates(
         'dither_spikes': dither_spikes,
         'jitter_spikes': jitter_spikes,
         'randomise_spikes': randomise_spikes,
-        'shuffle_isis': shuffle_isis}
+        'shuffle_isis': shuffle_isis,
+        'joint_isi_dithering': None}
 
     if surr_method not in surrogate_types.keys():
         raise ValueError('specified surr_method (=%s) not valid' % surr_method)
@@ -520,3 +975,5 @@ def surrogates(
     elif surr_method in ['randomise_spikes', 'shuffle_isis']:
         return surrogate_types[surr_method](
             spiketrain, n=n, decimals=decimals)
+    elif surr_method == 'joint_isi_dithering':
+        return JointISI(spiketrain).dithering(n)


Test patch:
diff --git a/elephant/test/test_conversion.py b/elephant/test/test_conversion.py
index 4d9db02..f0610b7 100644
--- a/elephant/test/test_conversion.py
+++ b/elephant/test/test_conversion.py
@@ -212,6 +212,22 @@ class BinnedSpikeTrainTestCase(unittest.TestCase):
         self.assertEqual(binned_matrix.get_num_of_spikes(),
                          sum(n_spikes_per_row_from_indices))
 
+    def test_get_num_of_spikes(self):
+        spiketrains = [self.spiketrain_a, self.spiketrain_b]
+        for spiketrain in spiketrains:
+            binned = cv.BinnedSpikeTrain(spiketrain, num_bins=10,
+                                         binsize=1 * pq.s, t_start=0 * pq.s)
+            self.assertEqual(binned.get_num_of_spikes(),
+                             len(binned.spike_indices[0]))
+        binned_matrix = cv.BinnedSpikeTrain(spiketrains, num_bins=10,
+                                            binsize=1 * pq.s)
+        n_spikes_per_row = binned_matrix.get_num_of_spikes(axis=1)
+        n_spikes_per_row_from_indices = list(map(len,
+                                                 binned_matrix.spike_indices))
+        assert_array_equal(n_spikes_per_row, n_spikes_per_row_from_indices)
+        self.assertEqual(binned_matrix.get_num_of_spikes(),
+                         sum(n_spikes_per_row_from_indices))
+
     def test_binned_spiketrain_sparse(self):
         a = neo.SpikeTrain([1.7, 1.8, 4.3] * pq.s, t_stop=10.0 * pq.s)
         b = neo.SpikeTrain([1.7, 1.8, 4.3] * pq.s, t_stop=10.0 * pq.s)
diff --git a/elephant/test/test_kernels.py b/elephant/test/test_kernels.py
index a9191e5..dcc1a3f 100644
--- a/elephant/test/test_kernels.py
+++ b/elephant/test/test_kernels.py
@@ -69,8 +69,9 @@ class kernel_TestCase(unittest.TestCase):
                        kernel_type in self.kernel_types]
         for kernel in kernel_list:
             b = kernel.boundary_enclosing_area_fraction(self.fraction).magnitude
-            restric_defdomain = \
-                np.linspace(-b, b, 2*b/kernel_resolution.magnitude) * sigma.units
+            n_points = int(2 * b / kernel_resolution.magnitude)
+            restric_defdomain = np.linspace(
+                -b, b, num=n_points) * sigma.units
             kern = kernel(restric_defdomain)
             norm = spint.cumtrapz(y=kern.magnitude,
                                   x=restric_defdomain.magnitude)[-1]
@@ -88,9 +89,9 @@ class kernel_TestCase(unittest.TestCase):
                            kernel_type in self.kernel_types]
             for kernel in kernel_list:
                 b = kernel.boundary_enclosing_area_fraction(self.fraction).magnitude
-                restric_defdomain = \
-                    np.linspace(-b, b, 2*b/kernel_resolution.magnitude) * \
-                    sigma.units
+                n_points = int(2 * b / kernel_resolution.magnitude)
+                restric_defdomain = np.linspace(
+                    -b, b, num=n_points) * sigma.units
                 kern = kernel(restric_defdomain)
                 av_integr = kern * restric_defdomain
                 average = spint.cumtrapz(y=av_integr.magnitude,
@@ -117,13 +118,14 @@ class kernel_TestCase(unittest.TestCase):
         for fraction in np.arange(0.15, 1.0, 0.4):
             for kernel in kernel_list:
                 b = kernel.boundary_enclosing_area_fraction(fraction).magnitude
-                restric_defdomain = \
-                    np.linspace(-b, b, 2*b/kernel_resolution.magnitude) * \
-                    sigma.units
+                n_points = int(2 * b / kernel_resolution.magnitude)
+                restric_defdomain = np.linspace(
+                    -b, b, num=n_points) * sigma.units
                 kern = kernel(restric_defdomain)
                 frac = spint.cumtrapz(y=kern.magnitude,
                                       x=restric_defdomain.magnitude)[-1]
                 self.assertAlmostEqual(frac, fraction, delta=0.002)
 
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/elephant/test/test_spike_train_surrogates.py b/elephant/test/test_spike_train_surrogates.py
index 046580e..7cf436a 100644
--- a/elephant/test/test_spike_train_surrogates.py
+++ b/elephant/test/test_spike_train_surrogates.py
@@ -8,15 +8,20 @@ unittests for spike_train_surrogates module.
 
 import unittest
 import elephant.spike_train_surrogates as surr
+import elephant.spike_train_generation as stg
 import numpy as np
+from numpy.testing import assert_array_almost_equal
 import quantities as pq
 import neo
-
-np.random.seed(0)
+import random
 
 
 class SurrogatesTestCase(unittest.TestCase):
 
+    def setUp(self):
+        np.random.seed(0)
+        random.seed(0)
+
     def test_dither_spikes_output_format(self):
 
         st = neo.SpikeTrain([90, 150, 180, 350] * pq.ms, t_stop=5 * pq.s)
@@ -54,7 +59,7 @@ class SurrogatesTestCase(unittest.TestCase):
 
         np.random.seed(42)
         dither_values = np.random.random_sample((nr_surr, len(st)))
-        expected_non_dithered = np.sum(dither_values==0)
+        expected_non_dithered = np.sum(dither_values == 0)
 
         observed_non_dithered = 0
         for surrog in surrs:
@@ -309,11 +314,85 @@ class SurrogatesTestCase(unittest.TestCase):
             self.assertEqual(len(surrog), len(st))
         self.assertTrue(len(surrs2) == nr_surr2)
 
+    def test_joint_isi_dithering_format(self):
+
+        rate = 100.*pq.Hz
+        t_stop = 1.*pq.s
+        st = stg.homogeneous_poisson_process(rate, t_stop=t_stop)
+        n_surr = 2
+        dither = 10 * pq.ms
+
+        # Test fast version
+        joint_isi_instance = surr.JointISI(st, dither=dither)
+        surrs = joint_isi_instance.dithering(n_surr=n_surr)
+
+        self.assertIsInstance(surrs, list)
+        self.assertEqual(len(surrs), n_surr)
+        self.assertEqual(joint_isi_instance.method, 'fast')
+
+        for surrog in surrs:
+            self.assertIsInstance(surrog, neo.SpikeTrain)
+            self.assertEqual(surrog.units, st.units)
+            self.assertEqual(surrog.t_start, st.t_start)
+            self.assertEqual(surrog.t_stop, st.t_stop)
+            self.assertEqual(len(surrog), len(st))
+
+        # Test window_version
+        joint_isi_instance = surr.JointISI(st,
+                                           method='window',
+                                           dither=2*dither,
+                                           num_bins=50)
+        surrs = joint_isi_instance.dithering(n_surr=n_surr)
+
+        self.assertIsInstance(surrs, list)
+        self.assertEqual(len(surrs), n_surr)
+        self.assertEqual(joint_isi_instance.method, 'window')
+
+        for surrog in surrs:
+            self.assertIsInstance(surrog, neo.SpikeTrain)
+            self.assertEqual(surrog.units, st.units)
+            self.assertEqual(surrog.t_start, st.t_start)
+            self.assertEqual(surrog.t_stop, st.t_stop)
+            self.assertEqual(len(surrog), len(st))
+
+        # Test surrogate methods wrapper
+        surrs = surr.surrogates(
+            st, n=n_surr, surr_method='joint_isi_dithering')
+        self.assertIsInstance(surrs, list)
+        self.assertEqual(len(surrs), n_surr)
+
+        for surrog in surrs:
+            self.assertIsInstance(surrog, neo.SpikeTrain)
+            self.assertEqual(surrog.units, st.units)
+            self.assertEqual(surrog.t_start, st.t_start)
+            self.assertEqual(surrog.t_stop, st.t_stop)
+            self.assertEqual(len(surrog), len(st))
+
+    def test_joint_isi_dithering_empty_train(self):
+        st = neo.SpikeTrain([] * pq.ms, t_stop=500 * pq.ms)
+        surrog = surr.JointISI(st).dithering()[0]
+        self.assertEqual(len(surrog), 0)
+
+    def test_joint_isi_dithering_output(self):
+        st = stg.homogeneous_poisson_process_with_refr_period(
+            rate=100. * pq.Hz,
+            refr_period=3 * pq.ms,
+            t_stop=0.1 * pq.s)
+        surrog_st = surr.JointISI(st).dithering()[0]
+        ground_truth = [27.99123087, 33.92680035, 39.94386642, 46.60932914,
+                        51.77647236, 58.00172645,
+                        68.15027025, 75.79592778, 84.09942906, 97.34738152]
+        assert_array_almost_equal(surrog_st.magnitude, ground_truth)
+
+
+
+
 
 def suite():
     suite = unittest.makeSuite(SurrogatesTestCase, 'test')
     return suite
 
+
 if __name__ == "__main__":
     runner = unittest.TextTestRunner(verbosity=2)
     runner.run(suite())