group_erp_tf_GO.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mne\n",
    "import glob\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "# mne.viz.set_browser_backend('qt')\n",
    "#%matplotlib qt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "project_path = \"/project/def-emayada/q1k/experimental/hsj/\"\n",
    "pylossless_path = \"derivatives/pylossless/\"\n",
    "postproc_path = \"derivatives/postproc/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditions = ['dtbc','dtgc','dtoc']\n",
    "roi = ['VREF']\n",
    "decim = 2\n",
    "freqs = np.arange(2, 50, 2)\n",
    "n_cycles = freqs / 2\n",
    "\n",
    "epoch_files = glob.glob(project_path + pylossless_path + postproc_path + 'epoch_fif_files/GO/*epo.fif')\n",
    "for item in epoch_files:\n",
    "   print(\"[\", item, \"] \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch_files = glob.glob(project_path + pylossless_path + postproc_path + 'epoch_fif_files/GO/*epo.fif')\n",
    "averaging_dict = {labels:[] for labels in conditions}\n",
    "\n",
    "for filepath in epoch_files:\n",
    "    print(filepath)\n",
    "    new_epoch = mne.read_epochs(filepath)\n",
    "    # the following is to select the first 50% of epochs\n",
    "    #nepochs = len(new_epoch)\n",
    "    #print(\"NTrials = \" + str(nepochs))\n",
    "    #trial_cap = round(nepochs/2)\n",
    "    #print(\"trial cap = \" + str(trial_cap))\n",
    "    #new_epoch = new_epoch[:trial_cap]\n",
    "    for condition in conditions:\n",
    "        power, itc = mne.time_frequency.tfr_morlet(new_epoch[condition].pick(roi), n_cycles=n_cycles, return_itc=True, freqs=freqs, decim=decim)\n",
    "        averaging_dict[condition].append((new_epoch[condition].average(), power, itc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def condition_summary(condition_label):\n",
    "    print('Working on: ', condition_label)\n",
    "    grand_average = mne.grand_average([item[0] for item in averaging_dict[condition_label]])\n",
    "    display(grand_average)\n",
    "    grand_average.plot()\n",
    "    times = np.arange(0, .6, .05)\n",
    "    fig = grand_average.plot_topomap(times=times, colorbar=True)\n",
    "    fig.suptitle(condition_label)\n",
    "\n",
    "for condition in conditions:\n",
    "    condition_summary(condition)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_dict = {'dtbc':'gray', 'dtoc':'blue', 'dtgc':'red'}\n",
    "linestyle_dict = {'dtbc':'-', 'dtoc':'-', 'dtgc':'-'}\n",
    "\n",
    "evokeds = {\n",
    "    'dtbc': [item[0] for item in averaging_dict['dtbc']],\n",
    "    'dtoc': [item[0] for item in averaging_dict['dtoc']],\n",
    "    'dtgc': [item[0] for item in averaging_dict['dtgc']],\n",
    "}\n",
    "\n",
    "mne.viz.plot_compare_evokeds(evokeds,\n",
    "                             combine='mean',\n",
    "                             legend='lower right',\n",
    "                             picks=roi, show_sensors='upper right',\n",
    "                             colors=color_dict,\n",
    "                             linestyles=linestyle_dict,\n",
    "                             title='baseline vs. overlap vs. gap ERPs'\n",
    "                            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# True for ERSP, False for ITC\n",
    "def do_power_plotting(ersp=True):\n",
    "    indexer = 1 if ersp else 2\n",
    "    cond1 = mne.grand_average([item[indexer] for item in averaging_dict['dtbc']])\n",
    "    cond2 = mne.grand_average([item[indexer] for item in averaging_dict['dtgc']])\n",
    "\n",
    "    epochs_power_1 = np.array([item[indexer].data for item in averaging_dict['dtbc']])[:, 0, :, :]\n",
    "    epochs_power_2 = np.array([item[indexer].data for item in averaging_dict['dtgc']])[:, 0, :, :]\n",
    "\n",
    "    #times = 1e3 * epochs_power_1.times  # change unit to ms\n",
    "    times = 1e3 * averaging_dict['dtbc'][0][1].times\n",
    "    fig1, (ax1t, ax1b) = plt.subplots(2, 1, figsize=(6, 4))\n",
    "    fig1.subplots_adjust(0.12, 0.08, 0.96, 0.94, 0.2, 0.43)\n",
    "\n",
    "    ax1t.imshow(\n",
    "        epochs_power_1.mean(axis=0),\n",
    "        extent=[times[0], times[-1], freqs[0], freqs[-1]],\n",
    "        aspect=\"auto\",\n",
    "        origin=\"lower\",\n",
    "        cmap=\"RdBu_r\",\n",
    "    )\n",
    "\n",
    "    ax1b.imshow(\n",
    "        epochs_power_2.mean(axis=0),\n",
    "        extent=[times[0], times[-1], freqs[0], freqs[-1]],\n",
    "        aspect=\"auto\",\n",
    "        origin=\"lower\",\n",
    "        cmap=\"RdBu_r\",\n",
    "    )\n",
    "\n",
    "    ax1t.set_ylabel(\"Frequency (Hz)\")\n",
    "    ax1t.set_title(f\"target Induced power\")\n",
    "    ax1b.set_title(f\"target Inter Trial Coherence\")\n",
    "    ax1b.set_xlabel(\"Time (ms)\")\n",
    "\n",
    "\n",
    "    F_obs, clusters, cluster_p_values, H0 = mne.stats.permutation_cluster_test(\n",
    "        [epochs_power_1, epochs_power_2],\n",
    "        out_type=\"mask\",\n",
    "        n_permutations=100,\n",
    "        threshold=6.0,\n",
    "        tail=0,\n",
    "    ) # returns F difference, sampled, zscore\n",
    "\n",
    "    # Grab any ERSP type timing window and use it to scale to ms\n",
    "    times = 1e3 * averaging_dict['dtbc'][0][1].times\n",
    "\n",
    "    evoked_power_contrast = epochs_power_1.mean(axis=0) - epochs_power_2.mean(axis=0)\n",
    "    signs = np.sign(evoked_power_contrast)\n",
    "\n",
    "    F_obs_plot = np.nan * np.ones_like(F_obs)\n",
    "    for c, p_val in zip(clusters, cluster_p_values):\n",
    "        if p_val <= 0.05:\n",
    "            F_obs_plot[c] = F_obs[c] * signs[c]\n",
    "    max_F = np.nanmax(abs(F_obs_plot))\n",
    "\n",
    "    fig, (ax, ax2) = plt.subplots(2, 1, figsize=(6, 4))\n",
    "    ax.imshow(\n",
    "        F_obs,\n",
    "        extent=[times[0], times[-1], freqs[0], freqs[-1]],\n",
    "        aspect=\"auto\",\n",
    "        origin=\"lower\",\n",
    "        cmap=\"gray\",\n",
    "    )\n",
    "\n",
    "    ax.imshow(\n",
    "        F_obs_plot,\n",
    "        extent=[times[0], times[-1], freqs[0], freqs[-1]],\n",
    "        aspect=\"auto\",\n",
    "        origin=\"lower\",\n",
    "        cmap=\"RdBu_r\",\n",
    "        vmin=-max_F,\n",
    "        vmax=max_F,\n",
    "    )\n",
    "    ax.set_xlabel(\"Time (ms)\")\n",
    "    ax.set_ylabel(\"Frequency (Hz)\")\n",
    "    ax.set_title(\"Induced power\")\n",
    "\n",
    "    evoked_contrast = mne.combine_evoked(\n",
    "        [cond1, cond2], weights=[1, -1]\n",
    "    )\n",
    "    evoked_contrast.plot(axes=ax2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "do_power_plotting(ersp=True)\n",
    "do_power_plotting(ersp=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!jupyter nbconvert --output {\"session_reports/group_half-trials_vp.html\"} --TagRemovePreprocessor.remove_all_outputs_tags='{\"exclude\"}' --no-input --to html group_erp_tf_vp.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (Q1K-REST)",
   "language": "python",
   "name": "q1k_env"
  },
  "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}