EEG_CN_Python_Handbook_AdvancedAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3cace2b9",
   "metadata": {},
   "source": [
    "# **第四章: 高级脑电数据分析**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eb4e115",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da5615ba",
   "metadata": {},
   "source": [
    "高级脑电数据分析的这一章分成了以下四个部分：\n",
    "\n",
    "- **第一节: 批处理读取与存储数据**\n",
    "- **第二节: 基于分类的脑电解码**  \n",
    "- **第三节: 表征相似性分析**\n",
    "- **第四节: 反向编码模型**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73737a19",
   "metadata": {},
   "source": [
    "## 下载并导入需要用到的Python包"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d0eb1b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "! pip install inverted_encoding\n",
    "\n",
    "import numpy as np\n",
    "import sys\n",
    "import os\n",
    "from six.moves import urllib\n",
    "import gdown\n",
    "import zipfile\n",
    "import scipy.io as sio\n",
    "from scipy.stats import ttest_ind\n",
    "import h5py\n",
    "from tqdm import tqdm\n",
    "from neurora.decoding import tbyt_decoding_kfold, ct_decoding_kfold\n",
    "from neurora.rsa_plot import plot_tbyt_decoding_acc, \\\n",
    "                             plot_ct_decoding_acc, \\\n",
    "                             plot_rdm, plot_tbytsim_withstats\n",
    "from neurora.rdm_cal import eegRDM\n",
    "from neurora.corr_cal_by_rdm import rdms_corr\n",
    "from neurora.stuff import smooth_1d\n",
    "from inverted_encoding import IEM, permutation, circ_diff\n",
    "from decimal import Decimal\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams[\"figure.dpi\"] = 300"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa86a97b",
   "metadata": {},
   "source": [
    "# 第一节 - 批处理读取与存储数据"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0868cb89",
   "metadata": {},
   "source": [
    "### **预处理过的Demo数据2：**  \n",
    "（你将会在第二节、第三节和第四节中使用）\n",
    "原数据集基于Bae&Luck2019年发表在Journal of Neuroscience上的文章《Dissociable Decoding of Spatial Attention and Working Memory from EEG Oscillations and Sustained Potentials》实验2的数据  \n",
    "这是一个视觉工作记忆任务，被试要求记忆一个水滴形状的朝向，刺激呈现200ms  \n",
    "经过1300毫秒的延迟后呈现一个随机朝向的水滴形状，被试需要转动鼠标使朝向尽可能和记忆朝向一致  \n",
    "刺激可能有16种不同的朝向以及16种不同的位置  \n",
    "这里仅提取前5个被试的数据，数据为做过预处理并分好段之后的ERP数据（预处理参数见论文原文）并带有每一个试次的朝向和位置标签"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "315059c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 下载Demo数据2\n",
    "\n",
    "data_dir = \"data/\"\n",
    "if not os.path.exists(data_dir):\n",
    "  os.makedirs(data_dir)\n",
    "\n",
    "# 从Google Drive下载\n",
    "url = \"https://drive.google.com/file/d/1POBi0dckB00AKCIvpXHZFI1Oqd-jWRoj/view?usp=sharing\"\n",
    "filename = \"demo_data2\"\n",
    "filepath = data_dir + filename + \".zip\"\n",
    "\n",
    "# 下载数据\n",
    "gdown.download(url=url, output=filepath, quiet=False, fuzzy=True)\n",
    "print(\"Download completes!\")\n",
    "# 解压数据\n",
    "with zipfile.ZipFile(filepath, 'r') as zip:\n",
    "    zip.extractall(data_dir)\n",
    "print(\"Unzip completes!\")\n",
    "\n",
    "# 也可以通过百度网盘下载\n",
    "# 链接:https://pan.baidu.com/s/1mZdU9rL8vGn8-kjlUkF6PA  密码:8qih\n",
    "# 下载后解压，并移动到data文件夹下"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2cf1466",
   "metadata": {},
   "source": [
    "在'data'文件夹下包含5个被试的脑电数据.mat文件，在'labels'文件夹下包含每个被试所有试次对应刺激项的朝向标签文件和位置标签文件\n",
    "以被试编号'201'为例，'ori_201.txt'文件包含其每一个试次的记忆项的朝向信息，'pos_201.txt'文件包含其每一个试次的记忆项的位置信息\n",
    "数据中，导联数为27，采样率为250Hz，从刺激呈现前1.5s到刺激呈现后1.5s，每个试次包含750个时间点"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ab2f9f5",
   "metadata": {},
   "source": [
    "### 批处理读取Demo数据2并存成.h5文件"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f157f828",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5个被试的编号\n",
    "sub_ids = ['201', '202', '203', '204', '205']\n",
    "\n",
    "# 对5个被试的数据进行遍历\n",
    "for i, sub in enumerate(sub_ids):\n",
    "    \n",
    "    # 每个被试ERP数据、记忆项朝向标签和记忆项位置标签的文件地址\n",
    "    subdata_path = 'data/demo_data2/data/ERP' + sub + '.mat'\n",
    "    suborilabels_path = 'data/demo_data2/labels/ori_' + sub + '.txt'\n",
    "    subposlabels_path = 'data/demo_data2/labels/pos_' + sub + '.txt'\n",
    "    \n",
    "    # 读取ERP数据\n",
    "    subdata = sio.loadmat(subdata_path)['filtData']\n",
    "    # 读取记忆项的朝向和位置的标签\n",
    "    # 在.txt文件里，第一列为具体的朝向/位置值，第二列为16个朝向/位置对应的标签值（0-15的整数表示）\n",
    "    suborilabels = np.loadtxt(suborilabels_path)[:, 1]\n",
    "    subposlabels = np.loadtxt(subposlabels_path)[:, 1]\n",
    "    \n",
    "    # 打印被试编号、ERP数据矩阵的shape及两个标签矩阵的shape\n",
    "    print('sub' + sub)\n",
    "    print(subdata.shape)\n",
    "    print(suborilabels.shape)\n",
    "    print(subposlabels.shape)\n",
    "    \n",
    "    # 将每个被试的ERP数据及朝向和位置的标签以‘data’、'orilabels'和‘poslabels’\n",
    "    # 作为Keys存在一个.h5文件里\n",
    "    f = h5py.File('data/demo_data2/sub' + sub + '.h5', 'w')\n",
    "    f.create_dataset('data', data=subdata)\n",
    "    f.create_dataset('orilabels', data=suborilabels)\n",
    "    f.create_dataset('poslabels', data=subposlabels)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03e6be8c",
   "metadata": {},
   "source": [
    "# 第二节 - 基于分类的脑电解码"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "364a867e",
   "metadata": {},
   "source": [
    "### Get EEG Data and Labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e81a082e",
   "metadata": {},
   "outputs": [],
   "source": [
    "subs = [\"201\", \"202\", \"203\", \"204\", \"205\"]\n",
    "\n",
    "# 初始化三个变量data、label_ori和label_pos\n",
    "# 分别用于后续存储脑电数据、朝向的标签和位置的标签\n",
    "data = np.zeros([5, 640, 27, 500])\n",
    "label_ori = np.zeros([5, 640])\n",
    "label_pos = np.zeros([5, 640])\n",
    "\n",
    "for i, sub in enumerate(subs):\n",
    "\n",
    "    # 读取单个被试的.h5文件\n",
    "    subfile = h5py.File('data/demo_data2/sub' + sub + '.h5', 'r')\n",
    "    \n",
    "    # 获取EEG数据\n",
    "    subdata = subfile['data']\n",
    "    \n",
    "    # subdata的shape为[640, 27, 750]\n",
    "    # 640 - 试次数； 27 - 导联数； 750 - 时间点数（从-1.5s到1.5s）\n",
    "\n",
    "    # 取-0.5s到1.5s的数据\n",
    "    subdata = subdata[:, :, 250:]\n",
    "    \n",
    "    # 读取朝向和位置标签\n",
    "    sublabel_ori = subfile['orilabels']\n",
    "    sublabel_pos = subfile['poslabels']\n",
    "    \n",
    "    data[i] = subdata\n",
    "    label_ori[i] = sublabel_ori\n",
    "    label_pos[i] = sublabel_pos"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74c542dc",
   "metadata": {},
   "source": [
    "### 逐时间点解码与结果可视化"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e7798e8",
   "metadata": {},
   "source": [
    "#### 对朝向信息进行逐时间点解码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea08e2e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 这里使用NeuroRA的decoding模块下的tbyt_decoding_kfold()函数进行逐时间点解码\n",
    "# 传入数据（shape为[n_subs, n_trials, n_channels, n_times]）、对应每个试次的标签\n",
    "# 以及一些其他参数：\n",
    "# n=16 - 分类类别数为16（共16中朝向）\n",
    "# navg=13 - 对13个试次进行一次平均\n",
    "# time_win=5 - 时间窗为5（即5个时间点的数据进行一次平均）\n",
    "# time_step=5 - 时间步长为5（即每间隔5个时间点进行一次解码）\n",
    "# nfolds=3 & nrepeats=10 - 进行10次3折的交叉验证\n",
    "# smooth=True - 对结果进行平滑处理\n",
    "accs_ori = tbyt_decoding_kfold(data, label_ori, n=16, navg=13, time_win=5, time_step=5, \n",
    "                               nfolds=3, nrepeats=10, smooth=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edaf5cd1",
   "metadata": {},
   "source": [
    "#### 对朝向信息解码结果的可视化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33f2983c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 这里使用NeuroRA的decoding模块下的plot_tbyt_decoding_acc()函数进行统计分析与可视化\n",
    "# 其要求输入解码的结果（shape为[n_subs, n_results_times]）和一些其他参数：\n",
    "# start_time=-0.5 - 解码结果的开始时间为-0.5s\n",
    "# end_time=1.5 - 解码结果的结束时间为1.5s、\n",
    "# time_interval=0.02 - 两解码时间间隔为0.02s\n",
    "#（由于时间窗为5，数据原采样率为250Hz，5*(1/250)=0.02s）、\n",
    "# chance=0.0625 - 随机水平的分类准确率为0.0625（16个朝向，即1/16）\n",
    "# p=0.05 - p值的阈值、\n",
    "# cbpt=True - 进行cluster-based permutation test、\n",
    "# stats_time=[0, 1.5] - 统计0-1.5s的解码结果、\n",
    "# xlim=[-0.5, 1.5] - x轴显示范围从-0.5到1.5\n",
    "# ylim=[0.05, 0.15] - y轴显示范围从0到1.5\n",
    "plot_tbyt_decoding_acc(accs_ori, start_time=-0.5, end_time=1.5, time_interval=0.02, \n",
    "                       chance=0.0625, p=0.05, cbpt=True, stats_time=[0, 1.5], \n",
    "                       xlim=[-0.5, 1.5], ylim=[0.05, 0.15])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "316aef13",
   "metadata": {},
   "source": [
    "#### 对位置信息进行逐时间点解码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6d4891a",
   "metadata": {},
   "outputs": [],
   "source": [
    "accs_pos = tbyt_decoding_kfold(data, label_pos, n=16, navg=13, time_win=5, time_step=5, \n",
    "                               nfolds=3, nrepeats=10, smooth=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d47bfe31",
   "metadata": {},
   "source": [
    "#### 对位置信息解码结果的可视化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "194a6769",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_tbyt_decoding_acc(accs_pos, start_time=-0.5, end_time=1.5, time_interval=0.02, \n",
    "                       chance=0.0625, p=0.05, cbpt=True, stats_time=[0, 1.5], \n",
    "                       xlim=[-0.5, 1.5], ylim=[0.05, 0.15])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0402d0a",
   "metadata": {},
   "source": [
    "### 跨时域脑电解码"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75f1bf21",
   "metadata": {},
   "source": [
    "#### 对朝向和位置信息进行跨时域脑电解码"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7c163e9",
   "metadata": {},
   "source": [
    "朝向解码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0318ff81",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 这里使用NeuroRA的decoding模块下的ct_decoding_kfold()函数进行跨时域解码\n",
    "accs_crosstime_ori = ct_decoding_kfold(data, label_ori, n=16, navg=13,\n",
    "                                       time_win=5, time_step=5, \n",
    "                                       nfolds=3, nrepeats=10, smooth=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd28ba3d",
   "metadata": {},
   "source": [
    "位置解码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01cb615e",
   "metadata": {},
   "outputs": [],
   "source": [
    "accs_crosstime_pos = ct_decoding_kfold(data, label_pos, n=16, navg=13,\n",
    "                                       time_win=5, time_step=5, \n",
    "                                       nfolds=3, nrepeats=10, smooth=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2127cfff",
   "metadata": {},
   "source": [
    "#### 对朝向和位置跨时域解码结果的可视化"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8357aefc",
   "metadata": {},
   "source": [
    "朝向解码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77f3f519",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 这里使用NeuroRA的decoding模块下的plot_ct_decoding_acc()函数进行统计分析与可视化\n",
    "plot_ct_decoding_acc(accs_crosstime_ori, start_timex=-0.5, end_timex=1.5, \n",
    "                     start_timey=-0.5, end_timey=1.5,\n",
    "                     time_intervalx=0.02, time_intervaly=0.02, \n",
    "                     chance=0.0625, p=0.05, cbpt=True,\n",
    "                     stats_timex=[0, 1.5], stats_timey=[0, 1.5], \n",
    "                     xlim=[-0.5, 1.5], ylim=[-0.5, 1.5], clim=[0.06, 0.075])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b836e8bb",
   "metadata": {},
   "source": [
    "位置解码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18775a22",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_ct_decoding_acc(accs_crosstime_pos, start_timex=-0.5, end_timex=1.5, \n",
    "                     start_timey=-0.5, end_timey=1.5,\n",
    "                     time_intervalx=0.02, time_intervaly=0.02, \n",
    "                     chance=0.0625, p=0.05, cbpt=True,\n",
    "                     stats_timex=[0, 1.5], stats_timey=[0, 1.5], \n",
    "                     xlim=[-0.5, 1.5], ylim=[-0.5, 1.5], clim=[0.06, 0.075])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8ff2b46",
   "metadata": {},
   "source": [
    "# 第三节 - 表征相似性分析"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebd3e964",
   "metadata": {},
   "source": [
    "### 计算脑电表征不相似性矩阵（Representational Dissimilarity Matrices，RDMs）"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe0dbc2a",
   "metadata": {},
   "source": [
    "根据16个朝向和16个位置，分别计算逐时间点的16×16的表征朝向信息的脑电RDM和表征位置信息的脑电RDM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ae3e835",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_subs = 5\n",
    "n_trials = 640\n",
    "# 首先对数据进行重新划分\n",
    "# 分别获取16个朝向条件下的脑电数据和16个位置条件下的脑电数据\n",
    "# 初始化两个变量data_ori和data_pos分别用于存储朝向和位置的数据\n",
    "data_ori = np.zeros([16, n_subs, 40, 27, 500])\n",
    "data_pos = np.zeros([16, n_subs, 40, 27, 500])\n",
    "for sub in range(n_subs):\n",
    "    index_ori = np.zeros([16], dtype=int)\n",
    "    index_pos = np.zeros([16], dtype=int)\n",
    "    for i in range(n_trials):\n",
    "        ori = int(label_ori[sub, i])\n",
    "        pos = int(label_pos[sub, i])\n",
    "        data_ori[ori, sub, index_ori[ori]] = data[sub, i]\n",
    "        index_ori[ori] = index_ori[ori] + 1\n",
    "        data_pos[pos, sub, index_pos[pos]] = data[sub, i]\n",
    "        index_pos[pos] = index_pos[pos] + 1\n",
    "        \n",
    "# 使用NeuroRA的rdm_cal模块下的eegRDM()函数计算脑电RDMs\n",
    "# 朝向RDMs\n",
    "RDMs_ori = eegRDM(data_ori, sub_opt=1, chl_opt=0, time_opt=1, time_win=5, time_step=5)\n",
    "# 位置RDMs\n",
    "RDMs_pos = eegRDM(data_pos, sub_opt=1, chl_opt=0, time_opt=1, time_win=5, time_step=5)\n",
    "# 返回的RDMs_ori和RDMs_pos的shape均为[n_subs即5, n_results_time即100, 16, 16]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef582da4",
   "metadata": {},
   "source": [
    "### 构建基于假设的表征不相似矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae983c3a",
   "metadata": {},
   "source": [
    "这里基于两角度越接近相似性约高、两角度越不接近相似性越低的假设  \n",
    "构建了一个编码模型RDM，此例中，对朝向和位置信息均适用"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e1bfefb",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_RDM = np.zeros([16, 16])\n",
    "for i in range(16):\n",
    "    for j in range(16):\n",
    "        diff = np.abs(i - j)\n",
    "        if diff <= 8:\n",
    "            model_RDM[i, j] = diff / 8\n",
    "        else:\n",
    "            model_RDM[i, j] = (16 - diff) / 8\n",
    "\n",
    "# 使用NeuroRA的rsa_plot模块下的plot_rdm()函数可视化RDM\n",
    "conditions = [\"0°\", \"22.5°\", \"45°\", \"67.5°\", \"90°\", \"112.5°\", \"135°\", \"157.5°\", \"180°\",\n",
    "              \"202.5°\", \"225°\", \"247.5°\", \"270°\", \"292.5°\", \"315°\", \"337.5°\"]\n",
    "plot_rdm(model_RDM, percentile=True, conditions=conditions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6c626cc",
   "metadata": {},
   "source": [
    "### 相似性分析（Representational Similarity Analysis, RSA）"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9725d58",
   "metadata": {},
   "source": [
    "使用基于假设的模型RDM分别与朝向表征的脑电RDMs和位置表征的脑电RDMs计算相似性  \n",
    "来时序上追踪大脑何时符合对朝向的编码、何时符合对位置的编码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee82c7a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 使用NeuroRA的corr_cal_by_rdm模块下的rdms_corr()函数\n",
    "# 计算基于假设的RDM和朝向的EEG RDMs的相似性\n",
    "similarities_ori = rdms_corr(model_RDM, RDMs_ori)\n",
    "# 计算基于假设的RDM和位置的EEG RDMs的相似性\n",
    "similarities_pos = rdms_corr(model_RDM, RDMs_pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aa6c7e4",
   "metadata": {},
   "source": [
    "基于RSA的大脑对朝向信息编码的结果可视化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb0c5a36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 使用NeuroRA的rsa_plot模块下的plot_tbytsim_withstats()函数\n",
    "# 用法与plot_tbyt_decoding_acc()函数类似\n",
    "plot_tbytsim_withstats(similarities_ori, start_time=-0.5, end_time=1.5, \n",
    "                       time_interval=0.02, smooth=True, p=0.05, cbpt=True, \n",
    "                       stats_time=[0, 1.5], xlim=[-0.5, 1.5], ylim=[-0.1, 0.8])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a52172f",
   "metadata": {},
   "source": [
    "基于RSA的大脑对位置信息编码的结果可视化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a98d057",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_tbytsim_withstats(similarities_pos, start_time=-0.5, end_time=1.5, \n",
    "                       time_interval=0.02, smooth=True, p=0.05, cbpt=True, \n",
    "                       stats_time=[0, 1.5], xlim=[-0.5, 1.5], ylim=[-0.1, 0.8])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c510d91",
   "metadata": {},
   "source": [
    "# 第四节 - 反向编码模型"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d7ed29f",
   "metadata": {},
   "source": [
    "这部分，我们参考预印本“*Scotti, P. S., Chen, J., & Golomb, J. D. (2021). An enhanced inverted encoding model for neural reconstructions. bioRxiv*”使用强化版反向编码模型（enhanced inverted encoding model, eIEM）"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0165b2f5",
   "metadata": {},
   "source": [
    "#### 使用eIEM解码朝向信息"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e661f8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_subs = 5\n",
    "n_ts = 100\n",
    "\n",
    "# 初始化mae_ori来存储IEM解码的MAE结果\n",
    "mae_ori = np.zeros([n_subs, n_ts])\n",
    "\n",
    "# 使用Inverted_Encoding的IEM函数计算逐被试、逐时间点的MAE\n",
    "for t in tqdm(range(n_ts)):\n",
    "    for sub in range(n_subs):\n",
    "        \n",
    "        # 降采样数据 - 即对每5个时间点平均\n",
    "        data_t_sub = np.average(data[sub, :, :, t*5:t*5+5], axis=2)\n",
    "        \n",
    "        # 获取IEM预测结果\n",
    "        predictions, _, _, _ = IEM(data_t_sub, label_ori[sub].astype(int),\n",
    "                                   stim_max=16, nfolds=5, is_circular=True)\n",
    "        \n",
    "        # 计算MAE\n",
    "        mae_ori[sub, t] = np.mean(np.abs(circ_diff(predictions, label_ori[sub].astype(int), 16)))\n",
    "\n",
    "# 使用Inverted_Encoding的permutation()函数获取MAE的Null分布\n",
    "null_mae_ori = permutation(label_ori[sub].astype(int), stim_max=16,num_perm=5000, is_circular=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01f57020",
   "metadata": {},
   "source": [
    "#### 定义绘制eIEM结果的函数 - plot_iem_results()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb262819",
   "metadata": {},
   "source": [
    "为了可视化，我们计算$\\Delta$MAE（null MAE的均值减MAE）作为最终的IEM结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5af9767",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_iem_results(mae_ori, null_mae, start_time=0, end_time=1, time_interval=0.01,\n",
    "                     stats_time=[0, 1], xlim=[0, 1], ylim=[-0.1, 0.8]):\n",
    "\n",
    "    \"\"\"\n",
    "    参数：\n",
    "        mae：shape为[n_subs, n_times]的矩阵，对应每个被试的MAE结果\n",
    "        null_mae：shape为[num_perm]的矩阵，对应Null分布\n",
    "        start_time：起始时间，默认为0\n",
    "        end_time：结束时间，默认为1\n",
    "        time_interval：两时间点之间的间隔，默认为0.01\n",
    "        xlim：X轴范围，默认[0, 1]\n",
    "        ylim：Y轴范围，默认[0.4, 0.8]\n",
    "    \"\"\"\n",
    "\n",
    "    if len(np.shape(mae_ori)) != 2:\n",
    "\n",
    "        return \"Invalid input!\"\n",
    "\n",
    "    n = len(np.shape(mae_ori))\n",
    "\n",
    "    yminlim = ylim[0]\n",
    "    ymaxlim = ylim[1]\n",
    "\n",
    "    nsubs, nts = np.shape(mae_ori)\n",
    "    tstep = float(Decimal((end_time - start_time) / nts).quantize(Decimal(str(time_interval))))\n",
    "\n",
    "    if tstep != time_interval:\n",
    "        return \"Invalid input!\"\n",
    "\n",
    "    delta1 = (stats_time[0] - start_time) / tstep - int((stats_time[0] - start_time) / tstep)\n",
    "    delta2 = (stats_time[1] - start_time) / tstep - int((stats_time[1] - start_time) / tstep)\n",
    "\n",
    "    mae_ori = smooth_1d(mae_ori)\n",
    "\n",
    "    avg = np.average(np.average(null_mae)-mae_ori, axis=0)\n",
    "    err = np.zeros([nts])\n",
    "\n",
    "    for t in range(nts):\n",
    "        err[t] = np.std(mae_ori[:, t], ddof=1)/np.sqrt(nsubs)\n",
    "\n",
    "    ps = np.zeros([nts])\n",
    "    for t in range(nts):\n",
    "        ps[t] = ttest_ind(mae_ori[:, t], null_mae, alternative='less')[1]\n",
    "        if ps[t] < 0.05:\n",
    "            ps[t] = 1\n",
    "        else:\n",
    "            ps[t] = 0\n",
    "\n",
    "    print('\\nSignificant time-windows:')\n",
    "    for t in range(nts):\n",
    "        if t == 0 and ps[t] == 1:\n",
    "            print(str(int(start_time * 1000)) + 'ms to ', end='')\n",
    "        if t > 0 and ps[t] == 1 and ps[t - 1] == 0:\n",
    "            print(str(int((start_time + t * tstep) * 1000)) + 'ms to ', end='')\n",
    "        if t < nts - 1 and ps[t] == 1 and ps[t + 1] == 0:\n",
    "            print(str(int((start_time + (t + 1) * tstep) * 1000)) + 'ms')\n",
    "        if t == nts - 1 and ps[t] == 1:\n",
    "            print(str(int(end_time * 1000)) + 'ms')\n",
    "\n",
    "    for t in range(nts):\n",
    "        if ps[t] == 1:\n",
    "            plt.plot(t*tstep+start_time+0.5*tstep, (ymaxlim-yminlim)*0.95+yminlim, 's',\n",
    "                     color='r', alpha=0.8, markersize=2)\n",
    "            xi = [t*tstep+start_time, t*tstep+tstep+start_time]\n",
    "            ymin = [0]\n",
    "            ymax = [avg[t]-err[t]]\n",
    "            plt.fill_between(xi, ymax, ymin, facecolor='r', alpha=0.1)\n",
    "\n",
    "    fig = plt.gcf()\n",
    "    fig.set_size_inches(6.4, 3.6)\n",
    "\n",
    "    ax = plt.gca()\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "    ax.spines[\"left\"].set_linewidth(3)\n",
    "    ax.spines[\"left\"].set_position((\"data\", 0))\n",
    "    ax.spines[\"bottom\"].set_linewidth(3)\n",
    "    ax.spines['bottom'].set_position(('data', 0))\n",
    "    x = np.arange(start_time+0.5*tstep, end_time+0.5*tstep, tstep)\n",
    "    plt.plot(x, avg, color='r', alpha=0.95)\n",
    "    plt.fill_between(x, avg + err, avg - err, facecolor='r', alpha=0.75)\n",
    "    plt.ylim(yminlim, ymaxlim)\n",
    "    plt.xlim(xlim[0], xlim[1])\n",
    "    plt.tick_params(labelsize=12)\n",
    "    plt.xlabel('Time (s)', fontsize=16)\n",
    "    plt.ylabel(r'$\\Delta$MAE', fontsize=16)\n",
    "    \n",
    "    plt.show()\n",
    "\n",
    "    return ps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f05998fd",
   "metadata": {},
   "source": [
    "绘制朝向信息的IEM结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76d7c556",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_iem_results(mae_ori, null_mae_ori, start_time=-0.5, end_time=1.5, time_interval=0.02,\n",
    "                 xlim=[-0.5, 1.5], ylim=[-0.1, 0.6])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1ef4e01",
   "metadata": {},
   "source": [
    "#### 使用eIEM解码位置信息"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "697b2347",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 初始化mae_pos来存储IEM解码的MAE结果\n",
    "mae_pos = np.zeros([n_subs, n_ts])\n",
    "\n",
    "# 使用Inverted_Encoding的IEM函数计算逐被试、逐时间点的MAE\n",
    "for t in tqdm(range(n_ts)):\n",
    "    for sub in range(n_subs):\n",
    "        \n",
    "        # 降采样数据 - 即对每5个时间点平均\n",
    "        data_t_sub = np.average(data[sub, :, :, t*5:t*5+5], axis=2)\n",
    "        \n",
    "        # 获取IEM预测结果\n",
    "        predictions, _, _, _ = IEM(data_t_sub, label_pos[sub].astype(int), stim_max=16, nfolds=5, is_circular=True)\n",
    "        \n",
    "        # 计算MAE\n",
    "        mae_pos[sub, t] = np.mean(np.abs(circ_diff(predictions, label_pos[sub].astype(int), 16)))\n",
    "\n",
    "# 使用Inverted_Encoding的permutation()函数获取MAE的Null分布\n",
    "null_mae_pos = permutation(label_pos[sub].astype(int), stim_max=16, num_perm=5000, is_circular=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba9825af",
   "metadata": {},
   "source": [
    "绘制位置信息的IEM结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b053699b",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot_iem_results(mae_pos, null_mae_pos, start_time=-0.5, end_time=1.5, time_interval=0.02,\n",
    "                 xlim=[-0.5, 1.5], ylim=[-0.1, 0.6])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}