MEEG_preprocessing_and_concatenate_epochs_(object_oriented_programming).ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "MEEG preprocessing and concatenate epochs (object oriented programming).ipynb",
      "provenance": [],
      "authorship_tag": "ABX9TyMz5OBTMvZR+55NeelsTFiG",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/nmningmei/preprocessing_pipelines/blob/master/MEEG_preprocessing_and_concatenate_epochs_(object_oriented_programming).ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QlBuH4Y3hgNL",
        "outputId": "040c6ebf-7d17-4ee3-a3b0-ab7545bef7ca"
      },
      "source": [
        "!pip install -U mne\r\n",
        "!pip install -U autoreject\r\n",
        "\r\n",
        "# Install the PyDrive wrapper & import libraries.\r\n",
        "# This only needs to be done once per notebook.\r\n",
        "!pip install -U -q PyDrive\r\n",
        "from pydrive.auth import GoogleAuth\r\n",
        "from pydrive.drive import GoogleDrive\r\n",
        "from google.colab import auth\r\n",
        "from oauth2client.client import GoogleCredentials\r\n",
        "\r\n"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Collecting mne\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/17/a8/7d8a10345082d4807907a268016b52dfa869b0c412cd84aa1d1de86e1e39/mne-0.22.0-py3-none-any.whl (6.9MB)\n",
            "\u001b[K     |████████████████████████████████| 6.9MB 6.1MB/s \n",
            "\u001b[?25hRequirement already satisfied, skipping upgrade: scipy>=0.17.1 in /usr/local/lib/python3.7/dist-packages (from mne) (1.4.1)\n",
            "Requirement already satisfied, skipping upgrade: numpy>=1.11.3 in /usr/local/lib/python3.7/dist-packages (from mne) (1.19.5)\n",
            "Installing collected packages: mne\n",
            "Successfully installed mne-0.22.0\n",
            "Collecting autoreject\n",
            "  Downloading https://files.pythonhosted.org/packages/00/00/1d93f88be662a1a65cae78261aad6856740a173c6947d47deaf91d70c47c/autoreject-0.2.2-py3-none-any.whl\n",
            "Installing collected packages: autoreject\n",
            "Successfully installed autoreject-0.2.2\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z7L7XkGJihpg"
      },
      "source": [
        "# Authenticate and create the PyDrive client.\r\n",
        "# This only needs to be done once per notebook.\r\n",
        "auth.authenticate_user()\r\n",
        "gauth = GoogleAuth()\r\n",
        "gauth.credentials = GoogleCredentials.get_application_default()\r\n",
        "drive = GoogleDrive(gauth)\r\n",
        "\r\n",
        "# Download a file based on its file ID.\r\n",
        "#\r\n",
        "# A file ID looks like: laggVyWshwcyP6kEI-y_W3P8D26sz\r\n",
        "eeg_ids = {1:'https://drive.google.com/open?id=1KRAtEBkuwqFloxxB-UgmpgrdOihHMFUo'.split('id=')[-1],\r\n",
        "           2:'https://drive.google.com/open?id=1w5bqX4KkFZMFOVZkxFl8svP6ASCfPPx_'.split('id=')[-1],\r\n",
        "           3:'https://drive.google.com/open?id=1jz2e8i_yi6vHqycf_3fe4w89UvARVZar'.split('id=')[-1],\r\n",
        "           #4:'https://drive.google.com/open?id=1xuOURk9L0nSjQ00JCg8W_19o4twQ_hbP'.split('id=')[-1],\r\n",
        "          }\r\n",
        "vmrk_ids = {1:'https://drive.google.com/open?id=1_1FfJdY4UPNe-GdYPCZBFwVUJQxrUuEJ'.split('id=')[-1],\r\n",
        "            2:'https://drive.google.com/open?id=1JZiQ5NmTAtn7Fj3urLmEJNYL9RGWzcIs'.split('id=')[-1],\r\n",
        "            3:'https://drive.google.com/open?id=1hB9Pyqt-Dx6EoyyAIgkCaeBZEMeQ2aqx'.split('id=')[-1],\r\n",
        "            #4:'https://drive.google.com/open?id=1RiBZ97jxVxI6GrpoVpo1nAUikPA9B1t4'.split('id=')[-1],\r\n",
        "           }\r\n",
        "vhdr_ids = {1:'https://drive.google.com/open?id=1uOkwHe-dySh5chqCyppKxrzwb0fV2x17'.split('id=')[-1],\r\n",
        "            2:'https://drive.google.com/open?id=17kgiU8HDIIKH2-YYeTTwXrAMSR2pivlh'.split('id=')[-1],\r\n",
        "            3:'https://drive.google.com/open?id=1D90dGVYZD6FKmh18bHQN-xgyLUtjsI6C'.split('id=')[-1],\r\n",
        "            #4:'https://drive.google.com/open?id=1IBVNlU0NtAdSZnbqWLQL8Ffvk8i-Gm3q'.split('id=')[-1],\r\n",
        "           }"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HYajSekWh3zq"
      },
      "source": [
        "for key,value in eeg_ids.items():\r\n",
        "    downloaded = drive.CreateFile({'id':value})\r\n",
        "    downloaded.GetContentFile(f'Patxi_sesion{key}.eeg')\r\n",
        "for key,value in vmrk_ids.items():\r\n",
        "    downloaded = drive.CreateFile({'id':value})\r\n",
        "    downloaded.GetContentFile(f'Patxi_sesion{key}.vmrk')\r\n",
        "for key,value in vhdr_ids.items():\r\n",
        "    downloaded = drive.CreateFile({'id':value})\r\n",
        "    downloaded.GetContentFile(f'Patxi_sesion{key}.vhdr')\r\n",
        "        "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dkUHev2Ph9HO"
      },
      "source": [
        "from glob import glob\r\n",
        "import mne\r\n",
        "import numpy as np\r\n",
        "#from autoreject import (AutoReject,get_rejection_threshold)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "MKY4wK-XiAhA",
        "outputId": "d5dce514-07ed-4a37-b479-2ba4ab416431"
      },
      "source": [
        "!git clone https://github.com/nmningmei/preprocessing_pipelines"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Cloning into 'preprocessing_pipelines'...\n",
            "remote: Enumerating objects: 3, done.\u001b[K\n",
            "remote: Counting objects: 100% (3/3), done.\u001b[K\n",
            "remote: Compressing objects: 100% (3/3), done.\u001b[K\n",
            "remote: Total 54 (delta 0), reused 0 (delta 0), pack-reused 51\u001b[K\n",
            "Unpacking objects: 100% (54/54), done.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lJEi3mOCiKyG"
      },
      "source": [
        "!cp preprocessing_pipelines/EEG_preprocessing.py EEG_preprocessing.py "
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8PgGlIHJiPwR"
      },
      "source": [
        "import EEG_preprocessing"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bU-h_yqEi12D"
      },
      "source": [
        "# The \"FASTER\" algorithm"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eO5FcoL3idkh"
      },
      "source": [
        "import re\r\n",
        "import mne\r\n",
        "import scipy\r\n",
        "\r\n",
        "import numpy as np\r\n",
        "\r\n",
        "from scipy.stats import kurtosis\r\n",
        "from numpy import nanmean\r\n",
        "from mne.utils import logger\r\n",
        "\r\n",
        "def str2int(x):\r\n",
        "    if type(x) is str:\r\n",
        "        return float(re.findall(r'\\d+',x)[0])\r\n",
        "    else:\r\n",
        "        return x\r\n",
        "def find_outliers(X, threshold=3.0, max_iter=2):\r\n",
        "    \"\"\"Find outliers based on iterated Z-scoring.\r\n",
        " \r\n",
        "    This procedure compares the absolute z-score against the threshold.\r\n",
        "    After excluding local outliers, the comparison is repeated until no\r\n",
        "    local outlier is present any more.\r\n",
        "    \r\n",
        "    ########ATTENTION ATTENTION ATTENTION#####\r\n",
        "    # This function if removed from MNE-python code base\r\n",
        "    Parameters\r\n",
        "    ----------\r\n",
        "    X : np.ndarray of float, shape (n_elemenets,)\r\n",
        "        The scores for which to find outliers.\r\n",
        "    threshold : float\r\n",
        "        The value above which a feature is classified as outlier.\r\n",
        "    max_iter : int\r\n",
        "        The maximum number of iterations.\r\n",
        " \r\n",
        "    Returns\r\n",
        "    -------\r\n",
        "    bad_idx : np.ndarray of int, shape (n_features)\r\n",
        "        The outlier indices.\r\n",
        "    \"\"\"\r\n",
        "    from scipy.stats import zscore\r\n",
        "    my_mask = np.zeros(len(X), dtype=np.bool)\r\n",
        "    for _ in range(max_iter):\r\n",
        "        X = np.ma.masked_array(X, my_mask)\r\n",
        "        this_z = np.abs(zscore(X))\r\n",
        "        local_bad = this_z > threshold\r\n",
        "        my_mask = np.max([my_mask, local_bad], 0)\r\n",
        "        if not np.any(local_bad):\r\n",
        "            break\r\n",
        " \r\n",
        "    bad_idx = np.where(my_mask)[0]\r\n",
        "    return bad_idx\r\n",
        "def hurst(x):\r\n",
        "    \"\"\"Estimate Hurst exponent on a timeseries.\r\n",
        "    The estimation is based on the second order discrete derivative.\r\n",
        "    Parameters\r\n",
        "    ----------\r\n",
        "    x : 1D numpy array\r\n",
        "        The timeseries to estimate the Hurst exponent for.\r\n",
        "    Returns\r\n",
        "    -------\r\n",
        "    h : float\r\n",
        "        The estimation of the Hurst exponent for the given timeseries.\r\n",
        "    \"\"\"\r\n",
        "    y = np.cumsum(np.diff(x, axis=1), axis=1)\r\n",
        "\r\n",
        "    b1 = [1, -2, 1]\r\n",
        "    b2 = [1,  0, -2, 0, 1]\r\n",
        "\r\n",
        "    # second order derivative\r\n",
        "    y1 = scipy.signal.lfilter(b1, 1, y, axis=1)\r\n",
        "    y1 = y1[:, len(b1) - 1:-1]  # first values contain filter artifacts\r\n",
        "\r\n",
        "    # wider second order derivative\r\n",
        "    y2 = scipy.signal.lfilter(b2, 1, y, axis=1)\r\n",
        "    y2 = y2[:, len(b2) - 1:-1]  # first values contain filter artifacts\r\n",
        "\r\n",
        "    s1 = np.mean(y1 ** 2, axis=1)\r\n",
        "    s2 = np.mean(y2 ** 2, axis=1)\r\n",
        "\r\n",
        "    return 0.5 * np.log2(s2 / s1)\r\n",
        "\r\n",
        "def _freqs_power(data, sfreq, freqs):\r\n",
        "    fs, ps = scipy.signal.welch(data, sfreq,\r\n",
        "                                nperseg=2 ** int(np.log2(10 * sfreq) + 1),\r\n",
        "                                noverlap=0,\r\n",
        "                                axis=-1)\r\n",
        "    return np.sum([ps[..., np.searchsorted(fs, f)] for f in freqs], axis=0)\r\n",
        "\r\n",
        "def faster_bad_channels(epochs, picks=None, thres=3, use_metrics=None):\r\n",
        "    \"\"\"Implements the first step of the FASTER algorithm.\r\n",
        "    \r\n",
        "    This function attempts to automatically mark bad EEG channels by performing\r\n",
        "    outlier detection. It operated on epoched data, to make sure only relevant\r\n",
        "    data is analyzed.\r\n",
        "    Parameters\r\n",
        "    ----------\r\n",
        "    epochs : Instance of Epochs\r\n",
        "        The epochs for which bad channels need to be marked\r\n",
        "    picks : list of int | None\r\n",
        "        Channels to operate on. Defaults to EEG channels.\r\n",
        "    thres : float\r\n",
        "        The threshold value, in standard deviations, to apply. A channel\r\n",
        "        crossing this threshold value is marked as bad. Defaults to 3.\r\n",
        "    use_metrics : list of str\r\n",
        "        List of metrics to use. Can be any combination of:\r\n",
        "            'variance', 'correlation', 'hurst', 'kurtosis', 'line_noise'\r\n",
        "        Defaults to all of them.\r\n",
        "    Returns\r\n",
        "    -------\r\n",
        "    bads : list of str\r\n",
        "        The names of the bad EEG channels.\r\n",
        "    \"\"\"\r\n",
        "    metrics = {\r\n",
        "        'variance':    lambda x: np.var(x, axis=1),\r\n",
        "        'correlation': lambda x: nanmean(\r\n",
        "                           np.ma.masked_array(\r\n",
        "                               np.corrcoef(x),\r\n",
        "                               np.identity(len(x), dtype=bool)\r\n",
        "                           ),\r\n",
        "                           axis=0),\r\n",
        "        'hurst':       lambda x: hurst(x),\r\n",
        "        'kurtosis':    lambda x: kurtosis(x, axis=1),\r\n",
        "        'line_noise':  lambda x: _freqs_power(x, epochs.info['sfreq'],\r\n",
        "                                              [50, 60]),\r\n",
        "    }\r\n",
        "\r\n",
        "    if picks is None:\r\n",
        "        picks = mne.pick_types(epochs.info, meg=False, eeg=True, exclude=[])\r\n",
        "    if use_metrics is None:\r\n",
        "        use_metrics = metrics.keys()\r\n",
        "\r\n",
        "    # Concatenate epochs in time\r\n",
        "    data = epochs.get_data()\r\n",
        "    data = data.transpose(1, 0, 2).reshape(data.shape[1], -1)\r\n",
        "    data = data[picks]\r\n",
        "\r\n",
        "    # Find bad channels\r\n",
        "    bads = []\r\n",
        "    for m in use_metrics:\r\n",
        "        s = metrics[m](data)\r\n",
        "        b = [epochs.ch_names[picks[i]] for i in find_outliers(s, thres)]\r\n",
        "        logger.info('Bad by %s:\\n\\t%s' % (m, b))\r\n",
        "        bads.append(b)\r\n",
        "\r\n",
        "    return np.unique(np.concatenate(bads)).tolist()\r\n",
        "\r\n",
        "def _deviation(data):\r\n",
        "    \"\"\"Computes the deviation from mean for each channel in a set of epochs.\r\n",
        "    This is not implemented as a lambda function, because the channel means\r\n",
        "    should be cached during the computation.\r\n",
        "    \r\n",
        "    Parameters\r\n",
        "    ----------\r\n",
        "    data : 3D numpy array\r\n",
        "        The epochs (#epochs x #channels x #samples).\r\n",
        "    Returns\r\n",
        "    -------\r\n",
        "    dev : 1D numpy array\r\n",
        "        For each epoch, the mean deviation of the channels.\r\n",
        "    \"\"\"\r\n",
        "    ch_mean = np.mean(data, axis=2)\r\n",
        "    return ch_mean - np.mean(ch_mean, axis=0)\r\n",
        "\r\n",
        "def faster_bad_epochs(epochs, picks=None, thres=3, use_metrics=None):\r\n",
        "    \"\"\"Implements the second step of the FASTER algorithm.\r\n",
        "    \r\n",
        "    This function attempts to automatically mark bad epochs by performing\r\n",
        "    outlier detection.\r\n",
        "    Parameters\r\n",
        "    ----------\r\n",
        "    epochs : Instance of Epochs\r\n",
        "        The epochs to analyze.\r\n",
        "    picks : list of int | None\r\n",
        "        Channels to operate on. Defaults to EEG channels.\r\n",
        "    thres : float\r\n",
        "        The threshold value, in standard deviations, to apply. An epoch\r\n",
        "        crossing this threshold value is marked as bad. Defaults to 3.\r\n",
        "    use_metrics : list of str\r\n",
        "        List of metrics to use. Can be any combination of:\r\n",
        "            'amplitude', 'variance', 'deviation'\r\n",
        "        Defaults to all of them.\r\n",
        "    Returns\r\n",
        "    -------\r\n",
        "    bads : list of int\r\n",
        "        The indices of the bad epochs.\r\n",
        "    \"\"\"\r\n",
        "\r\n",
        "    metrics = {\r\n",
        "        'amplitude': lambda x: np.mean(np.ptp(x, axis=2), axis=1),\r\n",
        "        'deviation': lambda x: np.mean(_deviation(x), axis=1),\r\n",
        "        'variance':  lambda x: np.mean(np.var(x, axis=2), axis=1),\r\n",
        "    }\r\n",
        "\r\n",
        "    if picks is None:\r\n",
        "        picks = mne.pick_types(epochs.info, meg=False, eeg=True,\r\n",
        "                               exclude='bads')\r\n",
        "    if use_metrics is None:\r\n",
        "        use_metrics = metrics.keys()\r\n",
        "\r\n",
        "    data = epochs.get_data()[:, picks, :]\r\n",
        "\r\n",
        "    bads = []\r\n",
        "    for m in use_metrics:\r\n",
        "        s = metrics[m](data)\r\n",
        "        b = find_outliers(s, thres)\r\n",
        "        logger.info('Bad by %s:\\n\\t%s' % (m, b))\r\n",
        "        bads.append(b)\r\n",
        "\r\n",
        "    return np.unique(np.concatenate(bads)).tolist()\r\n",
        "\r\n",
        "def _power_gradient(ica, source_data):\r\n",
        "    # Compute power spectrum\r\n",
        "    f, Ps = scipy.signal.welch(source_data, ica.info['sfreq'])\r\n",
        "\r\n",
        "    # Limit power spectrum to upper frequencies\r\n",
        "    Ps = Ps[:, np.searchsorted(f, 25):np.searchsorted(f, 45)]\r\n",
        "\r\n",
        "    # Compute mean gradients\r\n",
        "    return np.mean(np.diff(Ps), axis=1)\r\n",
        "\r\n",
        "\r\n",
        "def faster_bad_components(ica, epochs, thres=3, use_metrics=None):\r\n",
        "    \"\"\"Implements the third step of the FASTER algorithm.\r\n",
        "    \r\n",
        "    This function attempts to automatically mark bad ICA components by\r\n",
        "    performing outlier detection.\r\n",
        "    Parameters\r\n",
        "    ----------\r\n",
        "    ica : Instance of ICA\r\n",
        "        The ICA operator, already fitted to the supplied Epochs object.\r\n",
        "    epochs : Instance of Epochs\r\n",
        "        The untransformed epochs to analyze.\r\n",
        "    thres : float\r\n",
        "        The threshold value, in standard deviations, to apply. A component\r\n",
        "        crossing this threshold value is marked as bad. Defaults to 3.\r\n",
        "    use_metrics : list of str\r\n",
        "        List of metrics to use. Can be any combination of:\r\n",
        "            'eog_correlation', 'kurtosis', 'power_gradient', 'hurst',\r\n",
        "            'median_gradient'\r\n",
        "        Defaults to all of them.\r\n",
        "    Returns\r\n",
        "    -------\r\n",
        "    bads : list of int\r\n",
        "        The indices of the bad components.\r\n",
        "    See also\r\n",
        "    --------\r\n",
        "    ICA.find_bads_ecg\r\n",
        "    ICA.find_bads_eog\r\n",
        "    \"\"\"\r\n",
        "    source_data = ica.get_sources(epochs).get_data().transpose(1,0,2)\r\n",
        "    source_data = source_data.reshape(source_data.shape[0], -1)\r\n",
        "\r\n",
        "    metrics = {\r\n",
        "        'eog_correlation': lambda x: x.find_bads_eog(epochs)[1],\r\n",
        "        'kurtosis':        lambda x: kurtosis(\r\n",
        "                               np.dot(\r\n",
        "                                   x.mixing_matrix_.T,\r\n",
        "                                   x.pca_components_[:x.n_components_]),\r\n",
        "                               axis=1),\r\n",
        "        'power_gradient':  lambda x: _power_gradient(x, source_data),\r\n",
        "        'hurst':           lambda x: hurst(source_data),\r\n",
        "        'median_gradient': lambda x: np.median(np.abs(np.diff(source_data)),\r\n",
        "                                               axis=1),\r\n",
        "        'line_noise':  lambda x: _freqs_power(source_data,\r\n",
        "                                              epochs.info['sfreq'], [50, 60]),\r\n",
        "    }\r\n",
        "\r\n",
        "    if use_metrics is None:\r\n",
        "        use_metrics = metrics.keys()\r\n",
        "\r\n",
        "    bads = []\r\n",
        "    for m in use_metrics:\r\n",
        "        scores = np.atleast_2d(metrics[m](ica))\r\n",
        "        for s in scores:\r\n",
        "            b = find_outliers(s, thres)\r\n",
        "            logger.info('Bad by %s:\\n\\t%s' % (m, b))\r\n",
        "            bads.append(b)\r\n",
        "\r\n",
        "    return np.unique(np.concatenate(bads)).tolist()\r\n",
        "\r\n",
        "def faster_bad_channels_in_epochs(epochs, picks=None, thres=3, use_metrics=None):\r\n",
        "    \"\"\"Implements the fourth step of the FASTER algorithm.\r\n",
        "    \r\n",
        "    This function attempts to automatically mark bad channels in each epochs by\r\n",
        "    performing outlier detection.\r\n",
        "    Parameters\r\n",
        "    ----------\r\n",
        "    epochs : Instance of Epochs\r\n",
        "        The epochs to analyze.\r\n",
        "    picks : list of int | None\r\n",
        "        Channels to operate on. Defaults to EEG channels.\r\n",
        "    thres : float\r\n",
        "        The threshold value, in standard deviations, to apply. An epoch\r\n",
        "        crossing this threshold value is marked as bad. Defaults to 3.\r\n",
        "    use_metrics : list of str\r\n",
        "        List of metrics to use. Can be any combination of:\r\n",
        "            'amplitude', 'variance', 'deviation', 'median_gradient'\r\n",
        "        Defaults to all of them.\r\n",
        "    Returns\r\n",
        "    -------\r\n",
        "    bads : list of lists of int\r\n",
        "        For each epoch, the indices of the bad channels.\r\n",
        "    \"\"\"\r\n",
        "\r\n",
        "    metrics = {\r\n",
        "        'amplitude':       lambda x: np.ptp(x, axis=2),\r\n",
        "        'deviation':       lambda x: _deviation(x),\r\n",
        "        'variance':        lambda x: np.var(x, axis=2),\r\n",
        "        'median_gradient': lambda x: np.median(np.abs(np.diff(x)), axis=2),\r\n",
        "        'line_noise':      lambda x: _freqs_power(x, epochs.info['sfreq'],\r\n",
        "                                                  [50, 60]),\r\n",
        "    }\r\n",
        "\r\n",
        "    if picks is None:\r\n",
        "        picks = mne.pick_types(epochs.info, meg=False, eeg=True,\r\n",
        "                               exclude='bads')\r\n",
        "    if use_metrics is None:\r\n",
        "        use_metrics = metrics.keys()\r\n",
        "\r\n",
        "    \r\n",
        "    data = epochs.get_data()[:, picks, :]\r\n",
        "\r\n",
        "    bads = [[] for i in range(len(epochs))]\r\n",
        "    for m in use_metrics:\r\n",
        "        s_epochs = metrics[m](data)\r\n",
        "        for i, s in enumerate(s_epochs):\r\n",
        "            b = [epochs.ch_names[picks[j]] for j in find_outliers(s, thres)]\r\n",
        "            logger.info('Epoch %d, Bad by %s:\\n\\t%s' % (i, m, b))\r\n",
        "            bads[i].append(b)\r\n",
        "\r\n",
        "    for i, b in enumerate(bads):\r\n",
        "        if len(b) > 0:\r\n",
        "            bads[i] = np.unique(np.concatenate(b)).tolist()\r\n",
        "            \r\n",
        "    return bads\r\n",
        "\r\n",
        "def _check_bad_channels(_epochs,picks,func = faster_bad_channels):\r\n",
        "        bad_channels_list = func(_epochs,picks = picks)\r\n",
        "        for ch_name in bad_channels_list:\r\n",
        "            _epochs.info['bads'].append(ch_name)\r\n",
        "        return bad_channels_list,_epochs"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nP1lRvDzi5hi"
      },
      "source": [
        "# The EEG preprocessing pipeline in \"Class\""
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Bjtgbz4Pi4Ti"
      },
      "source": [
        "class PreprocessingPipeline(object):\r\n",
        "    \"\"\"\r\n",
        "    Parameters\r\n",
        "    ------\r\n",
        "    raw : mne.io.Raw, raw continue EEG recording\r\n",
        "    events : numpy.narray, n_events X 3, events matrix\r\n",
        "    event_id : dict, event ids that match the last colun of the events matrix\r\n",
        "    notch_filter : int or float, default = 50\r\n",
        "        band for notch filtering\r\n",
        "    perform_ICA : bool, default = True\r\n",
        "        whether plan to perform ICA preprocessing on the data\r\n",
        "    highpass : int or float, default = None\r\n",
        "        if not None, apply high pass filtering on the data\r\n",
        "    highpass_ICA : int or float, default = 1.\r\n",
        "        highpass a copy of the raw data at 1. Hz for ICA because ICA is sensitive\r\n",
        "        to low-frequency drifts. We want to remove them before fitting the ICA\r\n",
        "    lowpass : int or float, default = None\r\n",
        "        if not None, apply low pass filtering on the data\r\n",
        "    tmin : int or float,default = -0.5\r\n",
        "        the starting time point for defining the Epochs\r\n",
        "    tmax : int or float, default = 2.5\r\n",
        "        the ending time point for defining the Epochs\r\n",
        "    baseline : tuple of int or float, default = (-0.5,0)\r\n",
        "        the baseline for defining the Epochs\r\n",
        "    interpolate_bad_channels : bool, default = True\r\n",
        "        whether to perform bad channels interpolation\r\n",
        "    \r\n",
        "    Examples\r\n",
        "    -----\r\n",
        "    >>> # first we initialize the pipeline as a data container\r\n",
        "    >>> pipeline        = PreprocessingPipeline(raw,\r\n",
        "    >>>                                         events,\r\n",
        "    >>>                                         event_id,\r\n",
        "    >>>                                         )\r\n",
        "    >>> # now we preprocess the data step by step\r\n",
        "    >>> pipeline.re_refernce()\r\n",
        "    >>> pipeline.notch_filtering()\r\n",
        "    >>> pipeline.filtering()\r\n",
        "    >>> pipeline.epoching()\r\n",
        "    >>> pipeline.mark_bad_channels()\r\n",
        "    >>> pipeline.mark_bad_epochs()\r\n",
        "    >>> pipeline.mark_bad_channels_for_each_epoch()\r\n",
        "    >>> pipeline.fit_ica()\r\n",
        "    >>> pipeline.mark_bad_ica_components_by_FASTER()\r\n",
        "    >>> pipeline.detect_artifacts()\r\n",
        "    >>> pipeline.apply_ica()\r\n",
        "    >>> pipeline.final_step()\r\n",
        "    >>> clean_epochs = piplein.clean_epochs\r\n",
        "    >>> for event_name,val in event_id.items():\r\n",
        "    >>>     evoked = clean_epochs[event_name].average()\r\n",
        "    >>>     evoked.plot_joint(title = event_name)\r\n",
        "    \"\"\"\r\n",
        "    def __init__(self,\r\n",
        "                 raw,\r\n",
        "                 events,\r\n",
        "                 event_id,\r\n",
        "                 notch_filter               = 50,\r\n",
        "                 perform_ICA                = True,\r\n",
        "                 highpass                   = None,\r\n",
        "                 highpass_ICA               = 1.,\r\n",
        "                 lowpass                    = None,\r\n",
        "                 tmin                       = -0.5,\r\n",
        "                 tmax                       = 2.5,\r\n",
        "                 baseline                   = (-.5,0),\r\n",
        "                 interpolate_bad_channels   = True,\r\n",
        "                 ):\r\n",
        "        super(PreprocessingPipeline,).__init__()\r\n",
        "        np.random.seed(12345)\r\n",
        "        \r\n",
        "        self.raw                            = raw\r\n",
        "        self.events                         = events\r\n",
        "        self.event_id                       = event_id\r\n",
        "        \r\n",
        "        self.notch_filter                   = notch_filter\r\n",
        "        \r\n",
        "        self.highpass                       = highpass\r\n",
        "        self.lowpass                        = lowpass\r\n",
        "        self.highpass_ICA                   = highpass_ICA\r\n",
        "        \r\n",
        "        self.tmin                           = tmin\r\n",
        "        self.tmax                           = tmax\r\n",
        "        self.baseline                       = baseline\r\n",
        "        \r\n",
        "        self.perform_ICA                    = perform_ICA\r\n",
        "        \r\n",
        "        self.interpolate_bad_channels       = interpolate_bad_channels\r\n",
        "    \"\"\"\r\n",
        "    necessary step: re-reference - explicitly\r\n",
        "    \"\"\"\r\n",
        "    def re_refernce(self,):\r\n",
        "        self.raw_ref ,_  = mne.set_eeg_reference(self.raw,\r\n",
        "                                                 ref_channels     = 'average',\r\n",
        "                                                 projection       = True,)\r\n",
        "        self.raw_ref.apply_proj() # it might tell you it already has been re-referenced, but do it anyway\r\n",
        "    \r\n",
        "    \"\"\"\r\n",
        "    necessary step: notch filtering\r\n",
        "    \"\"\"\r\n",
        "    def notch_filtering(self,):\r\n",
        "        notch_filter = self.notch_filter\r\n",
        "        # everytime before filtering, explicitly pick the type of channels you want\r\n",
        "        # to perform the filters\r\n",
        "        picks = mne.pick_types(self.raw_ref.info,\r\n",
        "                               meg = False, # No MEG\r\n",
        "                               eeg = True,  # YES EEG\r\n",
        "                               eog = self.perform_ICA,  # depends on ICA\r\n",
        "                               )\r\n",
        "        # regardless the bandpass filtering later, we should always filter\r\n",
        "        # for wire artifacts and their oscillations\r\n",
        "        self.raw_ref.notch_filter(np.arange(notch_filter,241,notch_filter),\r\n",
        "                                  picks = picks)\r\n",
        "    \r\n",
        "    \"\"\"\r\n",
        "    optional step: highpass, lowpass, bandpass filtering\r\n",
        "    \"\"\"\r\n",
        "    def filtering(self,):\r\n",
        "        lowpass = self.lowpass\r\n",
        "        highpass = self.highpass\r\n",
        "        \r\n",
        "        if np.logical_and(highpass is not None,lowpass is not None):\r\n",
        "            self.raw_ref = self.raw_ref.filter(highpass,lowpass)\r\n",
        "        else:\r\n",
        "            if lowpass is not None:\r\n",
        "                self.raw_ref = self.raw_ref.filter(None,lowpass,)\r\n",
        "            elif highpass is not None:\r\n",
        "                self.raw_ref = self.raw_ref.filter(highpass,None)\r\n",
        "        \r\n",
        "        self.raw_ref_for_ICA = self.raw_ref.copy().filter(self.highpass_ICA,lowpass)\r\n",
        "    \"\"\"\r\n",
        "    necessary step: epoching the raw, not-filtered data\r\n",
        "    \"\"\"\r\n",
        "    def epoching(self,detrend = 1,preload = True,):\r\n",
        "        picks = mne.pick_types(self.raw_ref.info,\r\n",
        "                               meg = False, # No MEG\r\n",
        "                               eeg = True,  # YES EEG\r\n",
        "                               eog = self.perform_ICA,  # depends on ICA\r\n",
        "                               )\r\n",
        "        self.epochs      = mne.Epochs(self.raw_ref,\r\n",
        "                                      self.events,    # numpy array\r\n",
        "                                      self.event_id,  # dictionary\r\n",
        "                                      tmin        = self.tmin,\r\n",
        "                                      tmax        = self.tmax,\r\n",
        "                                      baseline    = self.baseline, # range of time for computing the mean references for each channel and subtract these values from all the time points per channel\r\n",
        "                                      picks       = picks,\r\n",
        "                                      detrend     = detrend, # linear detrend\r\n",
        "                                      preload     = preload # must be true if we want to do further processing\r\n",
        "                                      )\r\n",
        "        self.epochs_for_ICA = mne.Epochs(self.raw_ref_for_ICA,\r\n",
        "                                         self.events,    # numpy array\r\n",
        "                                         self.event_id,  # dictionary\r\n",
        "                                         tmin        = self.tmin,\r\n",
        "                                         tmax        = self.tmax,\r\n",
        "                                         baseline    = self.baseline, # range of time for computing the mean references for each channel and subtract these values from all the time points per channel\r\n",
        "                                         picks       = picks,\r\n",
        "                                         detrend     = detrend, # linear detrend\r\n",
        "                                         preload     = preload # must be true if we want to do further processing\r\n",
        "                                         )\r\n",
        "    \"\"\"\r\n",
        "    optional step: mark bad channels, interpolate them in necessary\r\n",
        "    \"\"\"\r\n",
        "    \r\n",
        "    \r\n",
        "    def mark_bad_channels(self,\r\n",
        "                          check_epochs = True,\r\n",
        "                          check_epochs_for_ICA = True,\r\n",
        "                          ):\r\n",
        "        picks = mne.pick_types(self.raw_ref.info,\r\n",
        "                               meg = False, # No MEG\r\n",
        "                               eeg = True,  # YES EEG\r\n",
        "                               eog = self.perform_ICA,  # depends on ICA\r\n",
        "                               )\r\n",
        "        if check_epochs:\r\n",
        "            self.bad_channels_list_epochs,self.epochs = _check_bad_channels(self.epochs,picks)\r\n",
        "            if self.interpolate_bad_channels:\r\n",
        "                self.epochs.interpolate_bads()\r\n",
        "        if check_epochs_for_ICA:\r\n",
        "            self.bad_channels_list_epochs_for_ICA, self.epochs_for_ICA = _check_bad_channels(self.epochs_for_ICA,picks)\r\n",
        "            if self.interpolate_bad_channels:\r\n",
        "                self.epochs_for_ICA.interpolate_bads()\r\n",
        "    \"\"\"\r\n",
        "    optional step: mark bad epochs\r\n",
        "    \"\"\"\r\n",
        "    def mark_bad_epochs(self,\r\n",
        "                        check_epochs = True,\r\n",
        "                        check_epochs_for_ICA = True,\r\n",
        "                        ):\r\n",
        "        picks = mne.pick_types(self.raw_ref.info,\r\n",
        "                               meg = False, # No MEG\r\n",
        "                               eeg = True,  # YES EEG\r\n",
        "                               eog = self.perform_ICA,  # depends on ICA\r\n",
        "                               )\r\n",
        "        if check_epochs:\r\n",
        "            self.bad_epochs_list_epochs = faster_bad_epochs(self.epochs,picks)\r\n",
        "        if check_epochs_for_ICA:\r\n",
        "            self.bad_epochs_list_epochs_for_ICA = faster_bad_epochs(self.epochs_for_ICA,picks)\r\n",
        "    \"\"\"\r\n",
        "    optional step: mark bad channels for each epoch\r\n",
        "    \"\"\"\r\n",
        "    def mark_bad_channels_for_each_epoch(self,\r\n",
        "                                         ):\r\n",
        "        picks = mne.pick_types(self.raw_ref.info,\r\n",
        "                               meg = False, # No MEG\r\n",
        "                               eeg = True,  # YES EEG\r\n",
        "                               eog = self.perform_ICA,  # depends on ICA\r\n",
        "                               )\r\n",
        "        self.bad_channels_list_epochs,self.epochs = _check_bad_channels(\r\n",
        "                self.epochs,picks,func = faster_bad_channels_in_epochs)\r\n",
        "        if self.interpolate_bad_channels:\r\n",
        "            self.epochs.interpolate_bads()\r\n",
        "        \r\n",
        "        \r\n",
        "    \"\"\"\r\n",
        "    half necessary half optional step: ICA fitting, ##### Not applying the ICA yet #####\r\n",
        "    \"\"\"\r\n",
        "    def fit_ica(self,\r\n",
        "                method                  = dict(noise_cov = 'empirical',\r\n",
        "                                               ica       = 'fastica',),\r\n",
        "                rank                    = None,\r\n",
        "                n_components            = .99,\r\n",
        "                n_pca_components        = .99,\r\n",
        "                max_pca_components      = None,\r\n",
        "                max_iter                = int(3e3),\r\n",
        "                verbose                 = 1,\r\n",
        "                ):\r\n",
        "        picks = mne.pick_types(self.epochs_for_ICA.info,\r\n",
        "                               meg = False, # No MEG\r\n",
        "                               eeg = True,  # YES EEG\r\n",
        "                               eog = False,\r\n",
        "                               )\r\n",
        "        # calculate the noise covariance of the epochs\r\n",
        "        noise_cov   = mne.compute_covariance(self.epochs_for_ICA,\r\n",
        "                                             tmin                   = self.baseline[0],\r\n",
        "                                             tmax                   = self.baseline[1],\r\n",
        "                                             method                 = method['noise_cov'],\r\n",
        "                                             rank                   = rank,\r\n",
        "                                             )\r\n",
        "        # define an ica function\r\n",
        "        ica         = mne.preprocessing.ICA(n_components            = n_components,\r\n",
        "                                            n_pca_components        = n_pca_components,\r\n",
        "                                            max_pca_components      = max_pca_components,\r\n",
        "                                            method                  = method['ica'],\r\n",
        "                                            max_iter                = max_iter,\r\n",
        "                                            noise_cov               = noise_cov,\r\n",
        "                                            random_state            = 12345,\r\n",
        "                                            )\r\n",
        "        # fit the ica\r\n",
        "        ica.fit(self.epochs_for_ICA,\r\n",
        "                picks   = picks,\r\n",
        "                start   = self.tmin,\r\n",
        "                stop    = self.tmax,\r\n",
        "                decim   = 3,\r\n",
        "                tstep   = 1., # Length of data chunks for artifact rejection in seconds. It only applies if inst is of type Raw.\r\n",
        "                verbose = verbose, # change to False if you want to omit the print information\r\n",
        "                )\r\n",
        "        self.ica = ica\r\n",
        "    \"\"\"\r\n",
        "    optional step: mark bad ica components using FASTER algorithm\r\n",
        "    \"\"\"\r\n",
        "    def mark_bad_ica_components_by_FASTER(self,):\r\n",
        "        self.bad_ica_list_by_FASTER = faster_bad_components(self.ica,self.epochs_for_ICA)\r\n",
        "        for idx in self.bad_ica_list_by_FASTER:\r\n",
        "            self.ica.exclude.append(idx)\r\n",
        "    \"\"\"\r\n",
        "    optional step: detect bad ica components using MNE-python function\r\n",
        "    \"\"\"\r\n",
        "    def detect_artifacts(self,\r\n",
        "                         eog_ch         = ['TP9','TP10','PO9','PO10'],\r\n",
        "                         eog_criterion  = 0.4, # arbitary choice\r\n",
        "                         skew_criterion = 1,   # arbitary choice\r\n",
        "                         kurt_criterion = 1,   # arbitary choice\r\n",
        "                         var_criterion  = 1,   # arbitary choice\r\n",
        "                         ):\r\n",
        "        # search for artificial ICAs automatically\r\n",
        "        # most of these hyperparameters were used in a unrelated published study\r\n",
        "        self.ica.detect_artifacts(self.epochs_for_ICA,\r\n",
        "                                  eog_ch         = eog_ch,\r\n",
        "                                  eog_criterion  = eog_criterion,\r\n",
        "                                  skew_criterion = skew_criterion,\r\n",
        "                                  kurt_criterion = kurt_criterion,\r\n",
        "                                  var_criterion  = var_criterion,\r\n",
        "                                  )\r\n",
        "    \"\"\"\r\n",
        "    half necessary half optional step: apply the fitted ica to the raw epochs\r\n",
        "    this step is necessary only if ICA is fit\r\n",
        "    \"\"\"\r\n",
        "    def apply_ica(self,):\r\n",
        "        self.epochs = self.ica.apply(self.epochs,\r\n",
        "                                     exclude = self.ica.exclude,\r\n",
        "                                     )\r\n",
        "    \"\"\"\r\n",
        "    final step: remove EOG channels from the list\r\n",
        "    \"\"\"\r\n",
        "    def final_step(self,):\r\n",
        "        self.clearn_epochs = self.epochs_ica.copy().pick_types(eeg = True, eog = False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fpFAfAHSjUEa"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}