Merge branch 'release/v0.5.2'

docs/source/conf.py

@@ -60,9 +60,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = '0.5.1'
+version = '0.5.2'
 # The full version, including alpha/beta/rc tags.
-release = '0.5.1'
+release = '0.5.2'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

docs/source/index.rst

@@ -55,6 +55,7 @@ Contents:
    methods
    pyfcn
    utils
+   stats
    visu
 
 

docs/source/methods.rst

@@ -29,7 +29,7 @@ The first digit refer to the PAC method :
 
 * 1 - Mean Vector Length (MVL - Canolty, 2006)
 * 2 - Kullback-Leibler Divergence (KLD - Tort, 2010)
-* 3 - Heigth-ratio (HR - Lakata, 2005)
+* 3 - Heigth-ratio (HR - Lakatos, 2005)
 * 4 - Normalized Direct PAC (ndPAC - Ozkürt, 2012)
 * 5 - Phase Synchrony (PS - Lakata, 2005)
 * 6 - Generalized Linear Model (GLM - Lakata, 2005) [IN PROGRESS]

docs/source/pyfcn.rst

@@ -9,6 +9,7 @@ The first thing to do is to define a `pacobj`. Then, from this objects, you can
 * :ref:`pacfilt` : if you already extracted your phase and amplitude, use this method to compute pac directly on it. 
 * :ref:`filtpac` : this is the all in one method. Starting from fresh data, this method will extract phase and amplitude and return the PAC estimation.
 * :ref:`pp` : compute the preferred-phase.
+* :ref:`erpac` : compute the Event-Related Phase-Amplitude Coupling (Voytek et al. 2013).
 
 .. _pacobj:
 
@@ -53,6 +54,25 @@ Use the following method to filter and compute the PAC directly.
 Preferred-phase
 ~~~~~~~~~~~~~~~
 
-Compute the preferred-phase (PP, see `this PP example <https://github.com/EtienneCmb/tensorpac/tree/master/examples/13_PreferredPhase.py>`_). 
+Compute the preferred-phase (PP, see `this PP example <https://github.com/EtienneCmb/tensorpac/tree/master/examples/12_PreferredPhase.py>`_). 
+
+.. figure::  picture/PP_planningVsExec.png
+   :align:   center
+
+   Example of preferred-phase for a motor task. This show how the amplitude is distributed according to phase values during a motor planning phase and a motor execution phase.  
 
 .. automethod:: tensorpac.pac.Pac.pp
+
+.. _erpac:
+
+Event-Related Phase-Amplitude Coupling
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Compute the ERPAC (see `this ERPAC example <https://github.com/EtienneCmb/tensorpac/tree/master/examples/13_EventRelatedPAC.py>`_). 
+
+.. figure::  picture/erpac.png
+   :align:   center
+
+   Example of Event-Related Phase-Amplitude Coupling (ERPAC, Voytek et al. 2013)
+
+.. automethod:: tensorpac.pac.Pac.erpac

docs/source/stats.rst

@@ -0,0 +1,15 @@
+.. _stats:
+
+Statistics
+==========
+
+Pearson correlation
+-------------------
+
+.. autofunction:: tensorpac.stats.pearson
+
+
+Circular correlation
+--------------------
+
+.. autofunction:: tensorpac.stats.circ_corrcc

docs/source/visu.rst

@@ -5,9 +5,12 @@ Visualization
 
 If matplotlib is install on your system, the main Pac class also include several plotting functions :
 
+* :ref:`pltpac` : standard plotting function for 2D arrays.
 * :ref:`pltcomo` : comodulogram plot, for finding the coupling location *PAC(phase, amplitude)*
 * :ref:`plttri` : triangular plot for finding where each band [Start, Stop] and the bandwidth *PAC(Starting frequency, Ending frequency)*
 * :ref:`pltpp` : polar plot for the preferred-phase representation *PP(phase, binned amplitude)*
+* :ref:`pltshow` : show the figure
+* :ref:`pltsave` : save the figure
 
 
 plotting method that rely on Matplotlib. This script is also present in the `example <https://github.com/EtienneCmb/tensorpac/blob/master/examples/9_PlayWithColors.py>`_ folder.
@@ -78,6 +81,13 @@ Code example
     plt.show()
 
 
+.. _pltpac:
+
+2D plot
+-------
+
+.. automethod:: tensorpac.pac.Pac.pacplot
+
 .. _pltcomo:
 
 Comodulogram
@@ -108,4 +118,18 @@ Polar plot
 .. figure::  picture/PP_planningVsExec.png
    :align:   center
 
-.. automethod:: tensorpac.pac.Pac.triplot
+.. automethod:: tensorpac.pac.Pac.triplot
+
+.. _pltshow:
+
+Show
+----
+
+.. automethod:: tensorpac.pac.Pac.show
+
+.. _pltsave:
+
+Save
+----
+
+.. automethod:: tensorpac.pac.Pac.savefig

dvp/erpac.py

@@ -7,31 +7,39 @@
 from brainpipe.feature import erpac
 
 ndatasets = 300
-x1, tvec = PacSignals(fpha=10, famp=100, ndatasets=ndatasets, noise=1, tmax=2)
-x2, tvec = PacSignals(fpha=20, famp=150, ndatasets=ndatasets, noise=1, tmax=2)
-x = np.concatenate((x1, x2), axis=1) + 5.
+x1, tvec = PacSignals(fpha=10, famp=100, ndatasets=ndatasets, noise=3, npts=1000, dpha=10, damp=10)
+x2 = np.random.rand(ndatasets, 700)
+x = np.concatenate((x1, x2), axis=1)
+T = np.arange(x.shape[1]) / 1024
 
 print('DATASET : ', x.shape)
 
-fpha = [9, 11]
+fpha = [8, 12]
 famp = [90, 110]
 
 t = time()
-p = Pac(idpac=(6, 0, 0), fpha=fpha, famp=famp, dcomplex='hilbert', width=12, nblocks=10)
-xpac, spac = p.fit(x, x, axis=1, nperm=10, traxis=0, njobs=-1)
+p = Pac(fpha=fpha, famp=(60, 150, 5, 1), dcomplex='wavelet', width=12)
+pha = p.filter(1024, x, axis=1, ftype='phase')
+amp = p.filter(1024, x, axis=1, ftype='amplitude')
+print(pha.shape, amp.shape)
+xpac, pval = p.erpac(pha, amp, traxis=1)
 print(xpac.shape)
 # 0/0
 print('Tensorpac : ', time()-t)
 
 # t = time()
-# e = erpac(1024, 2048, pha_f=fpha, amp_f=famp)
+# e = erpac(1024, 1500, pha_f=fpha, amp_f=famp)
 # bpac = np.squeeze(e.get(x.T, x.T, n_perm=10))[0]
 # print('Brainpipe : ', time()-t)
 # print(bpac.shape)
 
 # 0/0
-plt.subplot(1, 2, 1)
-plt.plot(np.squeeze(xpac))
+# plt.subplot(1, 2, 1)
+plt.pcolormesh(T, p.yvec, np.squeeze(xpac), cmap='Spectral_r')
+plt.title('ERPAC IS IN THE HOUSE')
+plt.xlabel('Time')
+plt.ylabel('Amplitude')
+# plt.plot(np.squeeze(xpac))
 plt.axis('tight')
 
 # plt.subplot(1, 2, 2)

dvp/ndpearson.py

@@ -34,11 +34,12 @@ def pear(x, y, axis=-1):
     cov = xy - n * mu_xy
     return cov / np.einsum('i...j, k...j->ik...', s_x, s_y)
 
-# p2 = generate_correlation_map(a, b)
+p2 = generate_correlation_map(a, b)
 p3 = pear(a, b)
+print(p3.shape)
 
 # print(p.shape, p2.shape, p3.shape)
-# print(np.allclose(p2, p3))
+print(np.allclose(p, p3))
 
 def circ_corrcc(alpha, x):
     n = alpha.shape[-1]

examples/03_ArtifialBandwidthSignals.py

@@ -15,5 +15,5 @@
         dcomplex='wavelet', width=6)
 pac = np.squeeze(p.filterfit(1024, data, data, axis=1, nperm=10)[0])
 
-p.comodulogram(pac.mean(-1), vmin=0.)
+p.comodulogram(pac.mean(-1), vmin=0., title=str(p))
 plt.show()

examples/05_ComparePacMethods.py

@@ -30,7 +30,7 @@
 for i, k in enumerate([1, 2, 3, 4, 5]):
     # Change the pac method :
     p.idpac = (k, 0, 0)
-    print('-> PAC using '+p.method)
+    print('-> PAC using '+str(p))
     # Compute only the PAC without filtering :
     xpac, _ = p.fit(phases, amplitudes, axis=2)
     # Plot :

examples/13_EventRelatedPAC.py

@@ -0,0 +1,57 @@
+"""Example of how to compute the ERPAC (Voytek et al 2013).
+
+Event-Related Phase-Amplitude Coupling (ERPAC) do not measure PAC across time
+cycle but instead, across trials (just as proposed JP. Lachaux with the
+PLV/PLS). Measuring across trials enable to have a real-time estimation of PAC.
+Warning, depending on your data, even with tensor calculation the ERPAC is
+significantly slower. Don't worry, take a coffee.
+
+In this example, we generate a signal that have a 10<->100hz coupling the first
+1000 points, then, the 700 following points are noise.
+"""
+import numpy as np
+from tensorpac import Pac, PacSignals
+
+# Generate a 10<->100hz coupling :
+ndatasets = 300
+npts = 1000
+sf = 1024.
+x1, tvec = PacSignals(fpha=10, famp=100, ndatasets=ndatasets, noise=2,
+                      npts=npts, dpha=10, damp=10, sf=sf)
+# Generate noise and concatenate the coupling and the noise :
+x2 = np.random.rand(ndatasets, 700)
+x = np.concatenate((x1, x2), axis=1)  # Shape : (ntrials, npts)
+time = np.arange(x.shape[1]) / sf
+
+# Define a PAC object :
+p = Pac(fpha=[9, 11], famp=(60, 140, 5, 1), dcomplex='wavelet', width=12)
+
+# Extract the phase and the amplitude :
+pha = p.filter(sf, x, axis=1, ftype='phase')  # Shape : (npha, ntrials, npts)
+amp = p.filter(sf, x, axis=1, ftype='amplitude')  # Shape : (namp, ntrials, npts)
+
+# Compute the ERPAC and use the traxis to specify that the trial axis is the
+# first one :
+erpac, pval = p.erpac(pha, amp, traxis=1)
+
+# Remove unused dimensions :
+erpac, pval = np.squeeze(erpac), np.squeeze(pval)
+
+# Plot without p-values :
+p.pacplot(erpac, time, p.yvec, xlabel='Time (second)', cmap='Spectral_r',
+          ylabel='Amplitude frequency', title=str(p), cblabel='ERPAC',
+          vmin=0., rmaxis=True)
+
+# Plot with every non-significant values masked in gray :
+# p.pacplot(erpac, time, p.yvec, xlabel='Time (second)', cmap='Spectral_r',
+#           ylabel='Amplitude frequency', title='ERPAC example', vmin=0.,
+#           vmax=1., pvalues=pval, bad='lightgray', plotas='contour',
+#           cblabel='ERPAC')
+
+# Plot with significiendy levels :
+# p.pacplot(erpac, time, p.yvec, xlabel='Time (second)', cmap='Spectral_r',
+#           ylabel='Amplitude frequency', title='ERPAC example', vmin=0.,
+#           vmax=1., pvalues=pval, levels=[1e-20, 1e-10, 1e-2, 0.05],
+#           levelcmap='inferno', plotas='contour', cblabel='ERPAC')
+# p.savefig('erpac.png', dpi=300)
+p.show()

setup.py

@@ -6,7 +6,7 @@
 from pip.req import parse_requirements
 from optparse import Option
 
-__version__ = "0.5.1"
+__version__ = "0.5.2"
 NAME = 'Tensorpac'
 AUTHOR = "Etienne Combrisson"
 MAINTAINER = "Etienne Combrisson"