From dd96e0a20d49ad2c8e0ebc1b8d7629bbf798f62f Mon Sep 17 00:00:00 2001 From: truikes Date: Wed, 16 Dec 2020 14:10:29 +0100 Subject: [PATCH 1/5] Draft version of 'pairwise phase consistency' --- elephant/phase_analysis.py | 62 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/elephant/phase_analysis.py b/elephant/phase_analysis.py index aa3a73345..46b995ee2 100644 --- a/elephant/phase_analysis.py +++ b/elephant/phase_analysis.py @@ -182,3 +182,65 @@ def spike_triggered_phase(hilbert_transform, spiketrains, interpolate): for i, entry in enumerate(result_times): result_times[i] = pq.Quantity(entry, units=entry[0].units).flatten() return result_phases, result_amps, result_times + + +def pairwise_phase_consistency(phases): + """ + The Pairwise Phase Consistency is an improved measure of phase consistency/phase locking value, + accounting for bias due to low trial counts. + + Published in Vinck et al., 2010 (https://www.sciencedirect.com/science/article/pii/S1053811910000959). + + Parameters + ---------- + phases : np.ndarray or list of np.ndarray + Spike-triggered phases (output from spike_triggered_phase). PPC is computed per array + + Returns + ------- + PPC : np.float or list of np.float + Pairwise Phase Consistency + + """ + + # Convert inputs to lists + if not isinstance(phases, list): + assert isinstance(phases, np.ndarray), 'Input should be an 1D np.array with phases' + phases = [phases] + + # Check if all elements are arrays + for p in phases: + assert isinstance(p, np.ndarray), 'Input should be an 1D np.array with phases or a list of those' + assert p.ndim == 1, 'Phase arrays should be 1D (use .flatten())' + + result_ppc = [] + + for p in phases: + n = p.shape[0] # nr of trials + + # Compute the distance between each pair of phases using dot product + # Optimize computation time using array multiplications instead of for loops + p_cos = np.cos(p) + p_cos_2d = np.empty((n, n), np.float32) # Note: don't think we need 64 precision + np.copyto(p_cos_2d, p_cos) + + p_sin = np.sin(p) + p_sin_2d = np.empty((n, n), np.float32) + np.copyto(p_sin_2d, p_sin) + + # By doing the elementwise multiplication of this matrix with its transpose, we get + # the distance between phases for all possible pairs of elements in p + temp_result = np.multiply(p_cos_2d, p_cos_2d.T) + np.multiply(p_sin_2d, p_sin_2d.T) + + # Now average over all elements in temp_results (the diagonal are 1 and should not be + # included) + di = np.diag_indices(n) + temp_result[di] = 0 + + ppc = np.sum(temp_result) / (n*n - n) + + result_ppc.append(ppc) + + return result_ppc + + From 9ad8c67b7f8234878f245cb36bbf75887db77c73 Mon Sep 17 00:00:00 2001 From: dizcza Date: Wed, 16 Dec 2020 15:40:00 +0100 Subject: [PATCH 2/5] bibtex citation, doc style --- doc/bib/elephant.bib | 11 ++++ doc/reference/phase_analysis.rst | 9 +++- elephant/phase_analysis.py | 87 +++++++++++++++++++------------- 3 files changed, 70 insertions(+), 37 deletions(-) diff --git a/doc/bib/elephant.bib b/doc/bib/elephant.bib index 353d8ef8d..577b50bc1 100644 --- a/doc/bib/elephant.bib +++ b/doc/bib/elephant.bib @@ -135,3 +135,14 @@ @article{Rostami17_3 pages = {3}, doi = {10.5281/zenodo.583814} } + +@article{Vinck2010_51, + title={The pairwise phase consistency: a bias-free measure of rhythmic neuronal synchronization}, + author={Vinck, Martin and van Wingerden, Marijn and Womelsdorf, Thilo and Fries, Pascal and Pennartz, Cyriel MA}, + journal={Neuroimage}, + volume={51}, + number={1}, + pages={112--122}, + year={2010}, + publisher={Elsevier} +} diff --git a/doc/reference/phase_analysis.rst b/doc/reference/phase_analysis.rst index 070f40675..61dbfc338 100644 --- a/doc/reference/phase_analysis.rst +++ b/doc/reference/phase_analysis.rst @@ -3,4 +3,11 @@ Spike-triggered LFP phase ========================= .. automodule:: elephant.phase_analysis - :members: + +References +---------- + +.. bibliography:: ../bib/elephant.bib + :labelprefix: ph + :keyprefix: phase- + :style: unsrt diff --git a/elephant/phase_analysis.py b/elephant/phase_analysis.py index 46b995ee2..0b69a48be 100644 --- a/elephant/phase_analysis.py +++ b/elephant/phase_analysis.py @@ -2,6 +2,12 @@ """ Methods for performing phase analysis. +.. autosummary:: + :toctree: toctree/phase_analysis + + spike_triggered_phase + pairwise_phase_consistency + :copyright: Copyright 2014-2018 by the Elephant team, see `doc/authors.rst`. :license: Modified BSD, see LICENSE.txt for details. """ @@ -12,7 +18,8 @@ import quantities as pq __all__ = [ - "spike_triggered_phase" + "spike_triggered_phase", + "pairwise_phase_consistency" ] @@ -185,62 +192,70 @@ def spike_triggered_phase(hilbert_transform, spiketrains, interpolate): def pairwise_phase_consistency(phases): - """ - The Pairwise Phase Consistency is an improved measure of phase consistency/phase locking value, - accounting for bias due to low trial counts. + r""" + The Pairwise Phase Consistency (PPC) :cite:`phase-Vinck2010_51` is an + improved measure of phase consistency/phase locking value, accounting for + bias due to low trial counts. + + The PPC is computed according to Eq. 14 and 15 of the cited paper. - Published in Vinck et al., 2010 (https://www.sciencedirect.com/science/article/pii/S1053811910000959). + .. math:: + \text{PPC} = \frac{2}{N(N-1)} \sum_{j=1}^{N-1} \sum_{k=j+1}^N + f(\theta_j, \theta_k) + + wherein the function :math:`f` computes the dot product between two unit + vectors and is defined by + + .. math:: + f(\phi, \omega) = \cos(\phi) \cos(\omega) + \sin(\phi) \sin(\omega) Parameters ---------- phases : np.ndarray or list of np.ndarray - Spike-triggered phases (output from spike_triggered_phase). PPC is computed per array + Spike-triggered phases (output from :func:`spike_triggered_phase`). + PPC is computed per array. Returns ------- - PPC : np.float or list of np.float + result_ppc : list of float Pairwise Phase Consistency """ - # Convert inputs to lists if not isinstance(phases, list): - assert isinstance(phases, np.ndarray), 'Input should be an 1D np.array with phases' phases = [phases] # Check if all elements are arrays - for p in phases: - assert isinstance(p, np.ndarray), 'Input should be an 1D np.array with phases or a list of those' - assert p.ndim == 1, 'Phase arrays should be 1D (use .flatten())' + if not isinstance(phases, (list, tuple)): + raise TypeError("Input must be a list of 1D numpy arrays with phases") + for phase_array in phases: + if not isinstance(phase_array, np.ndarray): + raise TypeError("Each entry of the input list must be an 1D " + "numpy array with phases") + if phase_array.ndim != 1: + raise ValueError("Phase arrays must be 1D (use .flatten())") result_ppc = [] - for p in phases: - n = p.shape[0] # nr of trials + for phase_array in phases: + n_trials = phase_array.shape[0] # Compute the distance between each pair of phases using dot product - # Optimize computation time using array multiplications instead of for loops - p_cos = np.cos(p) - p_cos_2d = np.empty((n, n), np.float32) # Note: don't think we need 64 precision - np.copyto(p_cos_2d, p_cos) - - p_sin = np.sin(p) - p_sin_2d = np.empty((n, n), np.float32) - np.copyto(p_sin_2d, p_sin) - - # By doing the elementwise multiplication of this matrix with its transpose, we get - # the distance between phases for all possible pairs of elements in p - temp_result = np.multiply(p_cos_2d, p_cos_2d.T) + np.multiply(p_sin_2d, p_sin_2d.T) - - # Now average over all elements in temp_results (the diagonal are 1 and should not be - # included) - di = np.diag_indices(n) - temp_result[di] = 0 - - ppc = np.sum(temp_result) / (n*n - n) - + # Optimize computation time using array multiplications instead of for + # loops + p_cos_2d = np.tile(np.cos(phase_array), reps=(n_trials, 1)) + p_sin_2d = np.tile(np.sin(phase_array), reps=(n_trials, 1)) + + # By doing the element-wise multiplication of this matrix with its + # transpose, we get the distance between phases for all possible pairs + # of elements in 'phase' + dot_prod = np.multiply(p_cos_2d, p_cos_2d.T) + \ + np.multiply(p_sin_2d, p_sin_2d.T) + + # Now average over all elements in temp_results (the diagonal are 1 + # and should not be included) + np.fill_diagonal(dot_prod, 0) + ppc = 2 / (n_trials * n_trials - n_trials) * np.sum(dot_prod) result_ppc.append(ppc) return result_ppc - - From 78e2ed15e66b7f27d947a5fb2e2df1e29edde84c Mon Sep 17 00:00:00 2001 From: truikes Date: Thu, 17 Dec 2020 12:04:18 +0100 Subject: [PATCH 3/5] Consider each entry in the `phases` input as a 'trial' from the same spiketrain --- doc/authors.rst | 3 +- doc/bib/elephant.bib | 15 +++++++ elephant/phase_analysis.py | 83 +++++++++++++++++++++++--------------- 3 files changed, 67 insertions(+), 34 deletions(-) diff --git a/doc/authors.rst b/doc/authors.rst index 48e81b1c0..42c4f6db6 100644 --- a/doc/authors.rst +++ b/doc/authors.rst @@ -51,6 +51,7 @@ Do you want to contribute to Elephant? Please refer to the * Regimantas Jurkus [1] * Philipp Steigerwald [12] * Manuel Ciba [12] +* Thijs Ruikes [13] 1. Institute of Neuroscience and Medicine (INM-6), Computational and Systems Neuroscience & Institute for Advanced Simulation (IAS-6), Theoretical Neuroscience, Jülich Research Centre and JARA, Jülich, Germany 2. Unité de Neurosciences, Information et Complexité, CNRS UPR 3293, Gif-sur-Yvette, France @@ -64,5 +65,5 @@ Do you want to contribute to Elephant? Please refer to the 10. Instituto de Neurobiología, Universidad Nacional Autónoma de México, Mexico City, Mexico 11. Case Western Reserve University (CWRU), Cleveland, OH, USA 12. BioMEMS Lab, TH Aschaffenburg University of applied sciences, Germany - +13. Cognitive and Systems Neuroscience (CSN) group, University of Amsterdam, Amsterdam, Netherlands If we've somehow missed you off the list we're very sorry - please let us know. diff --git a/doc/bib/elephant.bib b/doc/bib/elephant.bib index 577b50bc1..da2202841 100644 --- a/doc/bib/elephant.bib +++ b/doc/bib/elephant.bib @@ -146,3 +146,18 @@ @article{Vinck2010_51 year={2010}, publisher={Elsevier} } + +@article {PMID:22187161, + Title = {Improved measures of phase-coupling between spikes and the Local Field Potential}, + Author = {Vinck, Martin and Battaglia, Francesco Paolo and Womelsdorf, Thilo and Pennartz, Cyriel}, + DOI = {10.1007/s10827-011-0374-4}, + Number = {1}, + Volume = {33}, + Month = {August}, + Year = {2012}, + Journal = {Journal of computational neuroscience}, + ISSN = {0929-5313}, + Pages = {53—75}, + URL = {https://europepmc.org/articles/PMC3394239}, +} + diff --git a/elephant/phase_analysis.py b/elephant/phase_analysis.py index 0b69a48be..9c23a8d33 100644 --- a/elephant/phase_analysis.py +++ b/elephant/phase_analysis.py @@ -137,8 +137,8 @@ def spike_triggered_phase(hilbert_transform, spiketrains, interpolate): # Find index into signal for each spike ind_at_spike = ( - (spiketrain[sttimeind] - hilbert_transform[phase_i].t_start) / - hilbert_transform[phase_i].sampling_period). \ + (spiketrain[sttimeind] - hilbert_transform[phase_i].t_start) / + hilbert_transform[phase_i].sampling_period). \ simplified.magnitude.astype(int) # Append new list to the results for this spiketrain @@ -149,19 +149,19 @@ def spike_triggered_phase(hilbert_transform, spiketrains, interpolate): # Step through all spikes for spike_i, ind_at_spike_j in enumerate(ind_at_spike): - if interpolate and ind_at_spike_j+1 < len(times): + if interpolate and ind_at_spike_j + 1 < len(times): # Get relative spike occurrence between the two closest signal # sample points # if z->0 spike is more to the left sample # if z->1 more to the right sample - z = (spiketrain[sttimeind[spike_i]] - times[ind_at_spike_j]) /\ + z = (spiketrain[sttimeind[spike_i]] - times[ind_at_spike_j]) / \ hilbert_transform[phase_i].sampling_period # Save hilbert_transform (interpolate on circle) p1 = np.angle(hilbert_transform[phase_i][ind_at_spike_j]) p2 = np.angle(hilbert_transform[phase_i][ind_at_spike_j + 1]) interpolation = (1 - z) * np.exp(np.complex(0, p1)) \ - + z * np.exp(np.complex(0, p2)) + + z * np.exp(np.complex(0, p2)) p12 = np.angle([interpolation]) result_phases[spiketrain_i].append(p12) @@ -191,13 +191,19 @@ def spike_triggered_phase(hilbert_transform, spiketrains, interpolate): return result_phases, result_amps, result_times -def pairwise_phase_consistency(phases): +def pairwise_phase_consistency(phases, method='ppc0'): r""" - The Pairwise Phase Consistency (PPC) :cite:`phase-Vinck2010_51` is an + The Pairwise Phase Consistency (PPC0) :cite:`phase-Vinck2010_51` is an improved measure of phase consistency/phase locking value, accounting for bias due to low trial counts. - The PPC is computed according to Eq. 14 and 15 of the cited paper. + PPC0 is computed according to Eq. 14 and 15 of the cited paper. + + An improved version of the PPC (PPC1) :cite: `PMID:22187161` computes angular + difference ony between pairs of spikes within trials. + + PPC1 is not implemented yet + .. math:: \text{PPC} = \frac{2}{N(N-1)} \sum_{j=1}^{N-1} \sum_{k=j+1}^N @@ -213,7 +219,10 @@ def pairwise_phase_consistency(phases): ---------- phases : np.ndarray or list of np.ndarray Spike-triggered phases (output from :func:`spike_triggered_phase`). - PPC is computed per array. + If phases is a list of arrays, each array is considered a trial + + method : str + 'ppc0' - compute PPC between all pairs of spikes Returns ------- @@ -235,27 +244,35 @@ def pairwise_phase_consistency(phases): if phase_array.ndim != 1: raise ValueError("Phase arrays must be 1D (use .flatten())") - result_ppc = [] - - for phase_array in phases: - n_trials = phase_array.shape[0] - - # Compute the distance between each pair of phases using dot product - # Optimize computation time using array multiplications instead of for - # loops - p_cos_2d = np.tile(np.cos(phase_array), reps=(n_trials, 1)) - p_sin_2d = np.tile(np.sin(phase_array), reps=(n_trials, 1)) - - # By doing the element-wise multiplication of this matrix with its - # transpose, we get the distance between phases for all possible pairs - # of elements in 'phase' - dot_prod = np.multiply(p_cos_2d, p_cos_2d.T) + \ - np.multiply(p_sin_2d, p_sin_2d.T) - - # Now average over all elements in temp_results (the diagonal are 1 - # and should not be included) - np.fill_diagonal(dot_prod, 0) - ppc = 2 / (n_trials * n_trials - n_trials) * np.sum(dot_prod) - result_ppc.append(ppc) - - return result_ppc + if method not in ['ppc0']: + raise ValueError('For method choose out of: ["ppc0"]') + + phase_array = np.hstack(phases) + n_trials = phase_array.shape[0] # 'spikes' are 'trials' as in paper + + # Compute the distance between each pair of phases using dot product + # Optimize computation time using array multiplications instead of for + # loops + p_cos_2d = np.tile(np.cos(phase_array), reps=(n_trials, 1)) # TODO: optimize memory usage + p_sin_2d = np.tile(np.sin(phase_array), reps=(n_trials, 1)) + + # By doing the element-wise multiplication of this matrix with its + # transpose, we get the distance between phases for all possible pairs + # of elements in 'phase' + dot_prod = np.multiply(p_cos_2d, p_cos_2d.T) + \ + np.multiply(p_sin_2d, p_sin_2d.T) + + # Now average over all elements in temp_results (the diagonal are 1 + # and should not be included) + np.fill_diagonal(dot_prod, 0) + + if method == 'ppc0': + # Note: each pair i,j is computed twice in dot_prod. do not + # multiply by 2. n_trial * n_trials - n_trials = nr of filled elements + # in dot_prod + ppc = np.sum(dot_prod) / (n_trials * n_trials - n_trials) + return ppc + + elif method == 'ppc1': + # TODO: remove all indices from the same trial + return From 70ee3c8dcd3fdc7792815ff2417b95b5be01fc0f Mon Sep 17 00:00:00 2001 From: truikes Date: Mon, 21 Dec 2020 12:49:01 +0100 Subject: [PATCH 4/5] updated phase_analysis.py: -changed np.tile to np.broadcast_to -added dtype=np.float32 added tests in test_phase_analysis.py --- elephant/phase_analysis.py | 10 +-- elephant/test/test_phase_analysis.py | 126 +++++++++++++++++++++++++++ 2 files changed, 131 insertions(+), 5 deletions(-) diff --git a/elephant/phase_analysis.py b/elephant/phase_analysis.py index 9c23a8d33..aa8b21d39 100644 --- a/elephant/phase_analysis.py +++ b/elephant/phase_analysis.py @@ -253,14 +253,14 @@ def pairwise_phase_consistency(phases, method='ppc0'): # Compute the distance between each pair of phases using dot product # Optimize computation time using array multiplications instead of for # loops - p_cos_2d = np.tile(np.cos(phase_array), reps=(n_trials, 1)) # TODO: optimize memory usage - p_sin_2d = np.tile(np.sin(phase_array), reps=(n_trials, 1)) + p_cos_2d = np.broadcast_to(np.cos(phase_array), (n_trials, n_trials)) + p_sin_2d = np.broadcast_to(np.sin(phase_array), (n_trials, n_trials)) # By doing the element-wise multiplication of this matrix with its # transpose, we get the distance between phases for all possible pairs # of elements in 'phase' - dot_prod = np.multiply(p_cos_2d, p_cos_2d.T) + \ - np.multiply(p_sin_2d, p_sin_2d.T) + dot_prod = np.multiply(p_cos_2d, p_cos_2d.T, dtype=np.float32) + \ + np.multiply(p_sin_2d, p_sin_2d.T, dtype=np.float32) # TODO: agree on using this precision or not # Now average over all elements in temp_results (the diagonal are 1 # and should not be included) @@ -270,7 +270,7 @@ def pairwise_phase_consistency(phases, method='ppc0'): # Note: each pair i,j is computed twice in dot_prod. do not # multiply by 2. n_trial * n_trials - n_trials = nr of filled elements # in dot_prod - ppc = np.sum(dot_prod) / (n_trials * n_trials - n_trials) + ppc = np.sum(dot_prod) / (n_trials * n_trials - n_trials) # TODO: handle nan's return ppc elif method == 'ppc1': diff --git a/elephant/test/test_phase_analysis.py b/elephant/test/test_phase_analysis.py index 8385e3cc7..32d8d4c42 100644 --- a/elephant/test/test_phase_analysis.py +++ b/elephant/test/test_phase_analysis.py @@ -202,5 +202,131 @@ def test_regression_269(self): self.assertEqual(len(phases_noint[0]), 2) +class PairwisePhaseConsistencyTestCase(unittest.TestCase): + + @classmethod + def setUpClass(cls): # Note: using setUp makes the class call this + # function per test, while this way the function is called only + # 1 time per TestCase, slightly more efficient (0.5s tough) + + # Same setup as SpikeTriggerePhaseTestCase + tlen0 = 100 * pq.s + f0 = 20. * pq.Hz + fs0 = 1 * pq.ms + t0 = np.arange( + 0, tlen0.rescale(pq.s).magnitude, + fs0.rescale(pq.s).magnitude) * pq.s + cls.anasig0 = AnalogSignal( + np.sin(2 * np.pi * (f0 * t0).simplified.magnitude), + units=pq.mV, t_start=0 * pq.ms, sampling_period=fs0) + + # Spiketrain with perfect locking + cls.st_perfect = SpikeTrain( + np.arange(50, tlen0.rescale(pq.ms).magnitude - 50, 50) * pq.ms, + t_start=0 * pq.ms, t_stop=tlen0) + + # Spiketrain with inperfect locking + cls.st_inperfect = SpikeTrain( + [100., 100.1, 100.2, 100.3, 100.9, 101.] * pq.ms, + t_start=0 * pq.ms, t_stop=tlen0) + + # Generate 2 'bursting' spiketrains, both locking on sinus period, + # but with different strengths + n_spikes = 3 # n spikes per burst + burst_interval = (1 / f0.magnitude) * pq.s + burst_start_times = np.arange( + 0, + tlen0.rescale('ms').magnitude, + burst_interval.rescale('ms').magnitude + ) + + # Spiketrain with strong locking + burst_freq_strong = 200. * pq.Hz # strongly locking unit + burst_spike_interval = (1 / burst_freq_strong.magnitude) * pq.s + st_in_burst = np.arange( + 0, + burst_spike_interval.rescale('ms').magnitude * n_spikes, + burst_spike_interval.rescale('ms').magnitude + ) + st = [st_in_burst + t_offset for t_offset in burst_start_times] + st = np.hstack(st) * pq.ms + cls.st_bursting_strong = SpikeTrain(st, + t_start=0 * pq.ms, + t_stop=tlen0 + ) + + # Spiketrain with weak locking + burst_freq_weak = 100. * pq.Hz # weak locking unit + burst_spike_interval = (1 / burst_freq_weak.magnitude) * pq.s + st_in_burst = np.arange( + 0, + burst_spike_interval.rescale('ms').magnitude * n_spikes, + burst_spike_interval.rescale('ms').magnitude + ) + st = [st_in_burst + t_offset for t_offset in burst_start_times] + st = np.hstack(st) * pq.ms + cls.st_bursting_weak = SpikeTrain(st, + t_start=0 * pq.ms, + t_stop=tlen0 + ) + + def test_perfect_locking(self): + phases, _, _ = elephant.phase_analysis.spike_triggered_phase( + elephant.signal_processing.hilbert(self.anasig0), + self.st_perfect, + interpolate=True + ) + # Pass input as single array + ppc0 = elephant.phase_analysis.pairwise_phase_consistency( + phases[0], method='ppc0' + ) + self.assertEqual(ppc0, 1) + self.assertIsInstance(ppc0, float) + + # Pass input as list of arrays + n_phases = int(phases[0].shape[0] / 2) + phases_cut = [phases[0][i * 2:i * 2 + 2] for i in range(n_phases)] + ppc0 = elephant.phase_analysis.pairwise_phase_consistency( + phases_cut, method='ppc0' + ) + self.assertEqual(ppc0, 1) + self.assertIsInstance(ppc0, float) + + def test_inperfect_locking(self): + phases, _, _ = elephant.phase_analysis.spike_triggered_phase( + elephant.signal_processing.hilbert(self.anasig0), + self.st_inperfect, + interpolate=True + ) + # Pass input as single array + ppc0 = elephant.phase_analysis.pairwise_phase_consistency( + phases[0], method='ppc0' + ) + self.assertLess(ppc0, 1) + self.assertIsInstance(ppc0, float) + + def test_strong_vs_weak_locking(self): + phases_weak, _, _ = elephant.phase_analysis.spike_triggered_phase( + elephant.signal_processing.hilbert(self.anasig0), + self.st_bursting_weak, + interpolate=True + ) + # Pass input as single array + ppc0_weak = elephant.phase_analysis.pairwise_phase_consistency( + phases_weak[0], method='ppc0' + ) + phases_strong, _, _ = elephant.phase_analysis.spike_triggered_phase( + elephant.signal_processing.hilbert(self.anasig0), + self.st_bursting_strong, + interpolate=True + ) + # Pass input as single array + ppc0_strong = elephant.phase_analysis.pairwise_phase_consistency( + phases_strong[0], method='ppc0' + ) + + self.assertLess(ppc0_weak, ppc0_strong) + + if __name__ == '__main__': unittest.main() From ac51ba9e3cbff75188847a055557f5872e88dcac Mon Sep 17 00:00:00 2001 From: dizcza Date: Tue, 12 Jan 2021 16:38:32 +0100 Subject: [PATCH 5/5] fixed the citation, removed todos, pep8 --- doc/bib/elephant.bib | 2 +- elephant/phase_analysis.py | 14 ++++++-------- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/doc/bib/elephant.bib b/doc/bib/elephant.bib index da2202841..32ed60da2 100644 --- a/doc/bib/elephant.bib +++ b/doc/bib/elephant.bib @@ -147,7 +147,7 @@ @article{Vinck2010_51 publisher={Elsevier} } -@article {PMID:22187161, +@article {Vinck2012_33, Title = {Improved measures of phase-coupling between spikes and the Local Field Potential}, Author = {Vinck, Martin and Battaglia, Francesco Paolo and Womelsdorf, Thilo and Pennartz, Cyriel}, DOI = {10.1007/s10827-011-0374-4}, diff --git a/elephant/phase_analysis.py b/elephant/phase_analysis.py index aa8b21d39..c0546254e 100644 --- a/elephant/phase_analysis.py +++ b/elephant/phase_analysis.py @@ -199,8 +199,8 @@ def pairwise_phase_consistency(phases, method='ppc0'): PPC0 is computed according to Eq. 14 and 15 of the cited paper. - An improved version of the PPC (PPC1) :cite: `PMID:22187161` computes angular - difference ony between pairs of spikes within trials. + An improved version of the PPC (PPC1) :cite:`phase-Vinck2012_33` computes + angular difference ony between pairs of spikes within trials. PPC1 is not implemented yet @@ -230,13 +230,11 @@ def pairwise_phase_consistency(phases, method='ppc0'): Pairwise Phase Consistency """ - # Convert inputs to lists - if not isinstance(phases, list): + if isinstance(phases, np.ndarray): phases = [phases] - - # Check if all elements are arrays if not isinstance(phases, (list, tuple)): raise TypeError("Input must be a list of 1D numpy arrays with phases") + for phase_array in phases: if not isinstance(phase_array, np.ndarray): raise TypeError("Each entry of the input list must be an 1D " @@ -260,7 +258,7 @@ def pairwise_phase_consistency(phases, method='ppc0'): # transpose, we get the distance between phases for all possible pairs # of elements in 'phase' dot_prod = np.multiply(p_cos_2d, p_cos_2d.T, dtype=np.float32) + \ - np.multiply(p_sin_2d, p_sin_2d.T, dtype=np.float32) # TODO: agree on using this precision or not + np.multiply(p_sin_2d, p_sin_2d.T, dtype=np.float32) # Now average over all elements in temp_results (the diagonal are 1 # and should not be included) @@ -270,7 +268,7 @@ def pairwise_phase_consistency(phases, method='ppc0'): # Note: each pair i,j is computed twice in dot_prod. do not # multiply by 2. n_trial * n_trials - n_trials = nr of filled elements # in dot_prod - ppc = np.sum(dot_prod) / (n_trials * n_trials - n_trials) # TODO: handle nan's + ppc = np.sum(dot_prod) / (n_trials * n_trials - n_trials) return ppc elif method == 'ppc1':