Make install

Author: TomDLT

_downloads/plot_dar_model.py _downloads/05779cb76f55fa5acfb0044714058fd1/plot_dar_model.py

@@ -0,0 +1,122 @@
+"""
+Interface with MNE-python
+-------------------------
+This example shows how to use this package with MNE-python.
+
+It relies on the function ``raw_to_mask``, which takes as input a MNE.Raw
+instance and an events array, and returns the corresponding input signals
+and masks for the ``Comodulogram.fit`` method.
+"""
+import mne
+import numpy as np
+import os.path as op
+import matplotlib.pyplot as plt
+
+from pactools import simulate_pac, raw_to_mask, Comodulogram, MaskIterator
+
+###############################################################################
+# Simulate a dataset and save it out
+
+n_events = 100
+mu = 1.  # mean onset of PAC in seconds
+sigma = 0.25  # standard deviation of onset of PAC in seconds
+trial_len = 2.  # len of the simulated trial in seconds
+first_samp = 5  # seconds before the first sample and after the last
+
+fs = 200.  # Hz
+high_fq = 50.0  # Hz
+low_fq = 3.0  # Hz
+low_fq_width = 2.0  # Hz
+
+n_points = int(trial_len * fs)
+noise_level = 0.4
+
+
+def gaussian1d(array, mu, sigma):
+    return np.exp(-0.5 * ((array - mu) / sigma) ** 2)
+
+
+# leave one channel for events and make signal as long as events
+# with a bit of room on either side so events don't get cut off
+signal = np.zeros((2, int(n_points * n_events + 2 * first_samp * fs)))
+events = np.zeros((n_events, 3), dtype=int)
+events[:, 0] = np.arange((first_samp + mu) * fs,
+                         n_points * n_events + (first_samp + mu) * fs,
+                         trial_len * fs, dtype=int)
+events[:, 2] = np.ones((n_events))
+mod_fun = gaussian1d(np.arange(n_points), mu * fs, sigma * fs)
+for i in range(n_events):
+    signal_no_pac = simulate_pac(n_points=n_points, fs=fs,
+                                 high_fq=high_fq, low_fq=low_fq,
+                                 low_fq_width=low_fq_width,
+                                 noise_level=1.0, random_state=i)
+    signal_pac = simulate_pac(n_points=n_points, fs=fs,
+                              high_fq=high_fq, low_fq=low_fq,
+                              low_fq_width=low_fq_width,
+                              noise_level=noise_level, random_state=i)
+    signal[0, i * n_points:(i + 1) * n_points] = \
+        signal_pac * mod_fun + signal_no_pac * (1 - mod_fun)
+
+info = mne.create_info(['Ch1', 'STI 014'], fs, ['eeg', 'stim'])
+raw = mne.io.RawArray(signal, info)
+raw.add_events(events, stim_channel='STI 014')
+
+###############################################################################
+# Let's plot the signal and its power spectral density to visualize the data.
+# As shown in the plots below, there is a peak for the driver frequency at
+# 3 Hz and a peak for the carrier frequency at 50 Hz but phase-amplitude
+# coupling cannot be seen in the evoked plot by eye because the signal is
+# averaged over different phases for each epoch.
+
+raw.plot_psd(fmax=60)
+epochs = mne.Epochs(raw, events, tmin=-3, tmax=3)
+epochs.average().plot()
+
+###############################################################################
+# Let's save the raw object out for input/output demonstration purposes
+
+root = mne.utils._TempDir()
+raw.save(op.join(root, 'pac_example-raw.fif'))
+
+###############################################################################
+# Here we define how to build the epochs: which channels will be selected, and
+# on which time window around each event.
+
+raw = mne.io.Raw(op.join(root, 'pac_example-raw.fif'))
+events = mne.find_events(raw)
+
+# select the time interval around the events
+tmin, tmax = mu - 3 * sigma, mu + 3 * sigma
+# select the channels (phase_channel, amplitude_channel)
+ixs = (0, 0)
+
+###############################################################################
+# Then, we create the inputs with the function raw_to_mask, which creates the
+# input arrays and the mask arrays. These arrays are then given to a
+# comodulogram instance with the `fit` method, and the `plot` method draws the
+# results.
+
+# create the input array for Comodulogram.fit
+
+low_sig, high_sig, mask = raw_to_mask(raw, ixs=ixs, events=events, tmin=tmin,
+                                      tmax=tmax)
+
+###############################################################################
+# The mask is an iterable which goes over the _unique_ events in the event
+# array (if it is 3D). PAC is estimated where the `mask` is `False`.
+# Alternatively, we could also compute the `MaskIterator` object directly.
+# This is useful if you want to compute PAC on other kinds of time series,
+# for example source time courses.
+
+low_sig, high_sig = raw[ixs[0], :][0], raw[ixs[1], :][0]
+mask = MaskIterator(events, tmin, tmax, raw.n_times, raw.info['sfreq'])
+
+# create the instance of Comodulogram
+estimator = Comodulogram(fs=raw.info['sfreq'],
+                         low_fq_range=np.linspace(1, 10, 20), low_fq_width=2.,
+                         method='duprelatour', progress_bar=True)
+# compute the comodulogram
+estimator.fit(low_sig, high_sig, mask)
+# plot the results
+estimator.plot(tight_layout=False)
+plt.show()

_downloads/0789bcbe00d0db709615c2199c1448c7/plot_mne_epoch_masking.py

@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\nSurrogate analysis\n------------------\nThis example shows how to estimate a significance threshold in a comodulogram.\n\nA comodulogram shows the estimated PAC metric on a grid of frequency bands.\nIn absence of PAC, a PAC metric will return values close to zero, but not\nexactly zero. To estimate if a value is significantly far from zero, we use a\nsurrogate analysis.\n\nIn a surrogate analysis, we recompute the comodulogram many times, adding each\ntime a random time shift to remove any possible coupling between the\ncomponents. Nota that these time shifts have to be far from zero to effectively\nremove a potential coupling. These comodulograms give us an estimation of the\nfluctuation of the metric in the absence of PAC.\n\nTo derive a significance level from the list of comodulograms, we discuss here\ntwo different methods:\n- Computing a z-score on each couple of frequency, and thresholding at z = 4.\n- Computing a threshold at a given p-value, over a distribution of comodulogram\nmaxima.\n\n"
+        "\nSurrogate analysis\n------------------\nThis example shows how to estimate a significance threshold in a comodulogram.\n\nA comodulogram shows the estimated PAC metric on a grid of frequency bands.\nIn absence of PAC, a PAC metric will return values close to zero, but not\nexactly zero. To estimate if a value is significantly far from zero, we use a\nsurrogate analysis.\n\nIn a surrogate analysis, we recompute the comodulogram many times, adding each\ntime a random time shift to remove any possible coupling between the\ncomponents. Nota that these time shifts have to be far from zero to effectively\nremove a potential coupling. These comodulograms give us an estimation of the\nfluctuation of the metric in the absence of PAC.\n\nTo derive a significance level from the list of comodulograms, we discuss here\ntwo different methods:\n- Computing a z-score on each couple of frequency, and thresholding at z = 4.\n- Computing a threshold at a given p-value, over a distribution of comodulogram\nmaxima.\n"
       ]
     },
     {
@@ -161,7 +161,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.6.2"
+      "version": "3.8.1"
     }
   },
   "nbformat": 4,

_downloads/plot_surrogate_analysis.ipynb _downloads/097c1388cc7c939f7e77c63c1613cf11/plot_surrogate_analysis.ipynb

@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\nPhase shift and temporal delay in PAC\n-------------------------------------\n\nThis example disantangles the two distinct notions of phase shift and temporal\ndelay in phase-amplitude coupling. The phase shift ($\\phi$) is the phase of the\nslow oscillation which corresponds to the maximum amplitude of the fast\noscillation. The temporal delay ($\\tau$) is the delay between the two coupled\ncomponents. The two notions would be identical should the driver be a perfect\nstationary sinusoid.\n\n(1st line) When both are equal to zero, the high frequency bursts happen in\nthe driver's peaks.\n\n(2nd line) When $\\tau = 0$ and $\\phi != 0$, the bursts are shifted in time with\nrespect to the driver's peaks, and this shift varies depending on the\ninstantaneous frequency of the driver.\n\n(3rd line) When $\\tau != 0$ and $\\phi = 0$, the bursts are shifted in time with\nrespect to the driver's peaks, and this shift is constant over the signal. In\nthis case, note how the driver's phase corresponding to the bursts varies\ndepending on the instantaneous frequency of the driver.\n\n(4th line) Both $\\tau$ and $\\phi$ can also be both non-zero.\n\nThe temporal delay is estimated maximizing the likelihood on DAR models.\nThe phase shift is extracted from a DAR model fitted with the optimal\ntemporal delay.\n\nReferences\n==========\nDupre la Tour et al. (2017). Non-linear Auto-Regressive Models for\nCross-Frequency Coupling in Neural Time Series. bioRxiv, 159731.\n\n"
+        "\nPhase shift and temporal delay in PAC\n-------------------------------------\n\nThis example disantangles the two distinct notions of phase shift and temporal\ndelay in phase-amplitude coupling. The phase shift ($\\phi$) is the phase of the\nslow oscillation which corresponds to the maximum amplitude of the fast\noscillation. The temporal delay ($\\tau$) is the delay between the two coupled\ncomponents. The two notions would be identical should the driver be a perfect\nstationary sinusoid.\n\n(1st line) When both are equal to zero, the high frequency bursts happen in\nthe driver's peaks.\n\n(2nd line) When $\\tau = 0$ and $\\phi != 0$, the bursts are shifted in time with\nrespect to the driver's peaks, and this shift varies depending on the\ninstantaneous frequency of the driver.\n\n(3rd line) When $\\tau != 0$ and $\\phi = 0$, the bursts are shifted in time with\nrespect to the driver's peaks, and this shift is constant over the signal. In\nthis case, note how the driver's phase corresponding to the bursts varies\ndepending on the instantaneous frequency of the driver.\n\n(4th line) Both $\\tau$ and $\\phi$ can also be both non-zero.\n\nThe temporal delay is estimated maximizing the likelihood on DAR models.\nThe phase shift is extracted from a DAR model fitted with the optimal\ntemporal delay.\n\nReferences\n==========\nDupre la Tour et al. (2017). Non-linear Auto-Regressive Models for\nCross-Frequency Coupling in Neural Time Series. bioRxiv, 159731.\n"
       ]
     },
     {
@@ -46,7 +46,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.6.2"
+      "version": "3.8.1"
     }
   },
   "nbformat": 4,

_downloads/plot_phase_delay.ipynb _downloads/1698a7a11bdf113a93ec5f687e455c77/plot_phase_delay.ipynb

@@ -0,0 +1,162 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "%matplotlib inline"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\nInterface with MNE-python\n-------------------------\nThis example shows how to use this package with MNE-python.\n\nIt relies on the function ``raw_to_mask``, which takes as input a MNE.Raw\ninstance and an events array, and returns the corresponding input signals\nand masks for the ``Comodulogram.fit`` method.\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import mne\nimport numpy as np\nimport os.path as op\nimport matplotlib.pyplot as plt\n\nfrom pactools import simulate_pac, raw_to_mask, Comodulogram, MaskIterator"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Simulate a dataset and save it out\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "n_events = 100\nmu = 1.  # mean onset of PAC in seconds\nsigma = 0.25  # standard deviation of onset of PAC in seconds\ntrial_len = 2.  # len of the simulated trial in seconds\nfirst_samp = 5  # seconds before the first sample and after the last\n\nfs = 200.  # Hz\nhigh_fq = 50.0  # Hz\nlow_fq = 3.0  # Hz\nlow_fq_width = 2.0  # Hz\n\nn_points = int(trial_len * fs)\nnoise_level = 0.4\n\n\ndef gaussian1d(array, mu, sigma):\n    return np.exp(-0.5 * ((array - mu) / sigma) ** 2)\n\n\n# leave one channel for events and make signal as long as events\n# with a bit of room on either side so events don't get cut off\nsignal = np.zeros((2, int(n_points * n_events + 2 * first_samp * fs)))\nevents = np.zeros((n_events, 3), dtype=int)\nevents[:, 0] = np.arange((first_samp + mu) * fs,\n                         n_points * n_events + (first_samp + mu) * fs,\n                         trial_len * fs, dtype=int)\nevents[:, 2] = np.ones((n_events))\nmod_fun = gaussian1d(np.arange(n_points), mu * fs, sigma * fs)\nfor i in range(n_events):\n    signal_no_pac = simulate_pac(n_points=n_points, fs=fs,\n                                 high_fq=high_fq, low_fq=low_fq,\n                                 low_fq_width=low_fq_width,\n                                 noise_level=1.0, random_state=i)\n    signal_pac = simulate_pac(n_points=n_points, fs=fs,\n                              high_fq=high_fq, low_fq=low_fq,\n                              low_fq_width=low_fq_width,\n                              noise_level=noise_level, random_state=i)\n    signal[0, i * n_points:(i + 1) * n_points] = \\\n        signal_pac * mod_fun + signal_no_pac * (1 - mod_fun)\n\ninfo = mne.create_info(['Ch1', 'STI 014'], fs, ['eeg', 'stim'])\nraw = mne.io.RawArray(signal, info)\nraw.add_events(events, stim_channel='STI 014')"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Let's plot the signal and its power spectral density to visualize the data.\nAs shown in the plots below, there is a peak for the driver frequency at\n3 Hz and a peak for the carrier frequency at 50 Hz but phase-amplitude\ncoupling cannot be seen in the evoked plot by eye because the signal is\naveraged over different phases for each epoch.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "raw.plot_psd(fmax=60)\nepochs = mne.Epochs(raw, events, tmin=-3, tmax=3)\nepochs.average().plot()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Let's save the raw object out for input/output demonstration purposes\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "root = mne.utils._TempDir()\nraw.save(op.join(root, 'pac_example-raw.fif'))"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Here we define how to build the epochs: which channels will be selected, and\non which time window around each event.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "raw = mne.io.Raw(op.join(root, 'pac_example-raw.fif'))\nevents = mne.find_events(raw)\n\n# select the time interval around the events\ntmin, tmax = mu - 3 * sigma, mu + 3 * sigma\n# select the channels (phase_channel, amplitude_channel)\nixs = (0, 0)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Then, we create the inputs with the function raw_to_mask, which creates the\ninput arrays and the mask arrays. These arrays are then given to a\ncomodulogram instance with the `fit` method, and the `plot` method draws the\nresults.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# create the input array for Comodulogram.fit\n\nlow_sig, high_sig, mask = raw_to_mask(raw, ixs=ixs, events=events, tmin=tmin,\n                                      tmax=tmax)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The mask is an iterable which goes over the _unique_ events in the event\narray (if it is 3D). PAC is estimated where the `mask` is `False`.\nAlternatively, we could also compute the `MaskIterator` object directly.\nThis is useful if you want to compute PAC on other kinds of time series,\nfor example source time courses.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "low_sig, high_sig = raw[ixs[0], :][0], raw[ixs[1], :][0]\nmask = MaskIterator(events, tmin, tmax, raw.n_times, raw.info['sfreq'])\n\n# create the instance of Comodulogram\nestimator = Comodulogram(fs=raw.info['sfreq'],\n                         low_fq_range=np.linspace(1, 10, 20), low_fq_width=2.,\n                         method='duprelatour', progress_bar=True)\n# compute the comodulogram\nestimator.fit(low_sig, high_sig, mask)\n# plot the results\nestimator.plot(tight_layout=False)\nplt.show()"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.8.1"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

_downloads/1d03446eda919a1c40073436f5fe7156/plot_mne_epoch_masking.ipynb

@@ -79,7 +79,7 @@
 # Plug the model and the parameter grid into a GridSearchCV estimator
 # (GridSearchCVProgressBar is identical to GridSearchCV, but it adds a nice
 # progress bar to monitor progress.)
-gscv = GridSearchCVProgressBar(model, param_grid=param_grid,
+gscv = GridSearchCVProgressBar(model, param_grid=param_grid, cv=3,
                                return_train_score=False, verbose=1)
 
 # Fit the grid-search. We use `MultipleArray` to put together low_sig and