postprocess.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 导入必要的库"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import matlab.engine\n",
    "import time\n",
    "import scipy.io\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import scipy\n",
    "from accProcess import *\n",
    "sns.set_style('ticks')\n",
    "sns.set_context(\"poster\")\n",
    "plt.rcParams['font.sans-serif'] = 'Arial'\n",
    "# matlabeng = matlab.engine.start_matlab()   #启动matlab"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 指定模型和网络，加载数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_path = 'inp_arg3_model'\n",
    "result_path = 'basic_3inps_100class_logsa'\n",
    "figure_path = 'result_figs'\n",
    "fig_fmt = 'png'\n",
    "if not os.path.exists(os.path.join(model_path, result_path, figure_path)):\n",
    "    os.mkdir(os.path.join(model_path, result_path, figure_path))\n",
    "with open(os.path.join(model_path, 'data_preprocess_100class_logsa.pkl'), 'rb') as file:\n",
    "    dataset = pickle.loads(file.read())\n",
    "performance = scipy.io.loadmat(os.path.join(model_path, result_path, 'performance.mat'))\n",
    "results = scipy.io.loadmat(os.path.join(model_path, result_path, 'result.mat'))\n",
    "print(performance.keys())\n",
    "print(results.keys())\n",
    "period = np.logspace(np.log10(0.01), np.log10(10), 200)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train、Valid或Test数据集全部台站结果可视化"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 指定训练集，获取这些台站的数据，并计算残差"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取数据\n",
    "dataset_type = 'train'\n",
    "period_sec = [0.01, 0.1, 0.2, 0.5, 1.0, 2.0]\n",
    "data = []\n",
    "for i in range(len(dataset.train_data)):\n",
    "    data.append([])\n",
    "if dataset_type == 'train':\n",
    "    for i in range(len(dataset.train_data)):\n",
    "        data[i] = dataset.train_data[i]\n",
    "    label = dataset.train_label\n",
    "    pred = results['train_pred']\n",
    "elif dataset_type == 'valid':\n",
    "    for i in range(len(dataset.valid_data)):\n",
    "        data[i] = dataset.valid_data[i]\n",
    "    label = dataset.valid_label\n",
    "    pred = results['valid_pred']\n",
    "elif dataset_type == 'test':\n",
    "    for i in range(len(dataset.test_data)):\n",
    "        data[i] = dataset.test_data[i]\n",
    "    label = dataset.test_label\n",
    "    pred = results['test_pred']\n",
    "for i in range(len(data)):\n",
    "    data[i] = data[i] * dataset.data_std[i] + dataset.data_mean[i]\n",
    "label = label * dataset.label_std + dataset.label_mean\n",
    "\n",
    "# 计算残差\n",
    "resi_type = 'abs'\n",
    "if resi_type == 'abs':\n",
    "    residual = pred - label\n",
    "elif resi_type == 'rel':\n",
    "    residual = (pred - label) / label\n",
    "elif resi_type == 'log':\n",
    "    residual = np.log(pred) - np.log(label)\n",
    "\n",
    "data_scale = 'log'\n",
    "if data_scale == 'log':\n",
    "    label = np.log(10) * label\n",
    "    pred = np.log(10) * pred\n",
    "    data[0] = np.log(10) * data[0]\n",
    "    residual = np.log(10) * residual"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 不同周期下预测和实际反应谱值对比图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    r = np.corrcoef(label_points, pred_points)[0, 1]\n",
    "    xmin = np.min(label_points)\n",
    "    xmax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=label_points, y=pred_points, label='T=%.2fs' % T)\n",
    "    plt.plot([xmin, xmax], [xmin, xmax], 'k', linewidth=3)\n",
    "    plt.text(xmin ** 0.9 * xmax ** 0.1, xmin ** 0.1 * xmax ** 0.9, 'r=%.1f%%' % (100 * r), size=28)\n",
    "    plt.xlim([0.9 * xmin, 1.1 * xmax])\n",
    "    plt.ylim([0.9 * xmin, 1.1 * xmax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Target (gal)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Prediction (gal)')\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "    plt.legend(loc='lower right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_Sa_T_compare.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 绘制残差分布图"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随周期的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(period)\n",
    "num_bar = 10\n",
    "t_err, mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar), np.zeros(num_bar)\n",
    "for i in range(num_bar):\n",
    "    t_err[i] = period[int((i + 0.5) * N / num_bar)]\n",
    "    mean_err[i] = np.mean(residual[:, int(i * N / num_bar) : int((i + 1) * N / num_bar), 0])\n",
    "    std_err[i] = np.std(residual[:, int(i * N / num_bar) : int((i + 1) * N / num_bar), 0])\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "sns.scatterplot(x=np.tile(period, residual.shape[0]), y=residual[..., 0].ravel(), s=50, color='silver')\n",
    "plt.errorbar(t_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "plt.ylim([-0.6 * np.max(np.abs(residual[..., 0].ravel())), 0.6 * np.max(np.abs(residual[..., 0].ravel()))])\n",
    "plt.xscale('log')\n",
    "plt.xlabel('T (s)')\n",
    "plt.ylabel('Residual')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_residual_T.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随地表PGA的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载加速度时程信息数据\n",
    "event_msg = np.load('event_msg.npy', allow_pickle=True).item()\n",
    "# 获取数据集中每个地震动对应的PGA数据\n",
    "PGA_list = []\n",
    "for i in range(label.shape[0]):\n",
    "    if dataset_type == 'train':\n",
    "        name = dataset.train_events[dataset.train_idx[i]]\n",
    "    elif dataset_type == 'valid':\n",
    "        name = dataset.valid_events[dataset.valid_idx[i]]\n",
    "    elif dataset_type == 'test':\n",
    "        name = dataset.test_events[dataset.test_idx[i]]\n",
    "    PGA_list.append(981 * np.sqrt(event_msg[name][13] * event_msg[name][39]))\n",
    "\n",
    "value_list = np.array(PGA_list)\n",
    "\n",
    "# 将所有地震的PGA按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "num_bar = 10\n",
    "value_step = (np.log(value_max) - np.log(value_min)) / num_bar\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = np.exp((np.arange(num_bar) + 0.5) * value_step + np.log(value_min))\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if np.log(value) >= np.log(value_min) + i * value_step and np.log(value) < np.log(value_min) + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随PGA的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err[:-1], mean_err[:-1], std_err[:-1], fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('PGA (gal)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='lower right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_residual_PGA.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随$v_{s30}$的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载vs30数据并获取数据集中每条地震动对应的Vs30值\n",
    "vs30 = np.load('vs30.npy', allow_pickle=True).item()\n",
    "vs30_list = []\n",
    "for i in range(label.shape[0]):\n",
    "    if dataset_type == 'train':\n",
    "        name = dataset.train_events[dataset.train_idx[i]]\n",
    "    elif dataset_type == 'valid':\n",
    "        name = dataset.valid_events[dataset.valid_idx[i]]\n",
    "    elif dataset_type == 'test':\n",
    "        name = dataset.test_events[dataset.test_idx[i]]\n",
    "    vs30_list.append(vs30[name[:6]])\n",
    "value_list = np.array(vs30_list)\n",
    "\n",
    "# 将所有地震对应的Vs30按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.min(value_list)\n",
    "value_max = 1500\n",
    "num_bar = 10\n",
    "value_step = (value_max - value_min) / num_bar\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = value_min + (np.arange(num_bar) + 0.5) * value_step\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if value >= value_min + i * value_step and value < value_min + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随Vs30的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('$V_{s30}$ (m/s)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.legend(loc='upper right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_residual_vs30.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随震级的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震级列表\n",
    "value_list = data[2][:, 0, 0]\n",
    "\n",
    "# 将所有地震对应的震级按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.round(np.min(value_list)) - 0.5\n",
    "value_max = np.round(np.max(value_list)) + 0.5\n",
    "value_step = 1\n",
    "num_bar = int((value_max - value_min) / value_step)\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = np.arange(int(value_min + 0.5), int(value_max + 0.5))\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if value >= value_min + i * value_step and value < value_min + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随震级的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([value_min, value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Mg.')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.legend(loc='lower right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_residual_Mg.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随震源深度的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震源深度列表\n",
    "depth_list = data[2][:, 1, 0]\n",
    "depth_list[depth_list < 1.0] = 1.0\n",
    "value_list = depth_list\n",
    "\n",
    "# 将所有地震对应的震源深度按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "num_bar = 10\n",
    "value_step = (np.log(value_max) - np.log(value_min)) / num_bar\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = np.exp(np.log(value_min) + (np.arange(num_bar) + 0.5) * value_step)\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if np.log(value) >= np.log(value_min) + i * value_step and np.log(value) < np.log(value_min) + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随震源深度的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Depth (km)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='lower right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_residual_Depth.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随震中距的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震中距列表\n",
    "value_list = data[2][:, 2, 0]\n",
    "\n",
    "# 将所有地震对应的震中距按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "num_bar = 10\n",
    "value_step = (np.log(value_max) - np.log(value_min)) / num_bar\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = np.exp(np.log(value_min) + (np.arange(num_bar) + 0.5) * value_step)\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if np.log(value) >= np.log(value_min) + i * value_step and np.log(value) < np.log(value_min) + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随震中距的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('R (km)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='upper left')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_residual_Rup.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 绘制label和prediction的散点对比图"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随周期的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ymin = 0\n",
    "ymax = np.max([np.max(label[..., 0].ravel()), np.max(pred[..., 0].ravel())])\n",
    "plt.figure(figsize=(8, 6))\n",
    "sns.scatterplot(x=np.tile(period, label.shape[0]), y=label[..., 0].ravel(), s=20, label='Target')\n",
    "sns.scatterplot(x=np.tile(period, label.shape[0]), y=pred[..., 0].ravel(), s=20, label='Prediction')\n",
    "# plt.ylim([ymin, 0.2 * ymax])\n",
    "plt.xscale('log')\n",
    "plt.xlabel('T (s)')\n",
    "plt.ylabel('Sa (gal)')\n",
    "plt.legend()\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_Sa_scatter_T.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随地表PGA的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载加速度时程信息数据\n",
    "event_msg = np.load('event_msg.npy', allow_pickle=True).item()\n",
    "# 获取数据集中每个地震动对应的PGA数据\n",
    "PGA_list = []\n",
    "for i in range(label.shape[0]):\n",
    "    if dataset_type == 'train':\n",
    "        name = dataset.train_events[dataset.train_idx[i]]\n",
    "    elif dataset_type == 'valid':\n",
    "        name = dataset.valid_events[dataset.valid_idx[i]]\n",
    "    elif dataset_type == 'test':\n",
    "        name = dataset.test_events[dataset.test_idx[i]]\n",
    "    PGA_list.append(981 * np.sqrt(event_msg[name][13] * event_msg[name][39]))\n",
    "value_list = np.array(PGA_list)\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(value_min, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('PGA (gal)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='upper left')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_Sa_scatter_PGA.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随$v_{s30}$的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载vs30数据并获取数据集中每条地震动对应的Vs30值\n",
    "vs30 = np.load('vs30.npy', allow_pickle=True).item()\n",
    "vs30_list = []\n",
    "for i in range(label.shape[0]):\n",
    "    if dataset_type == 'train':\n",
    "        name = dataset.train_events[dataset.train_idx[i]]\n",
    "    elif dataset_type == 'valid':\n",
    "        name = dataset.valid_events[dataset.valid_idx[i]]\n",
    "    elif dataset_type == 'test':\n",
    "        name = dataset.test_events[dataset.test_idx[i]]\n",
    "    vs30_list.append(vs30[name[:6]])\n",
    "value_list = np.array(vs30_list)\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(0.8 * value_max, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('$V_{s30}$ (m/s)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.legend(loc='upper right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_Sa_scatter_vs30.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随震级的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震级列表\n",
    "value_list = data[2][:, 0, 0]\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(value_min, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Mg.')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.legend(loc='upper left')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_Sa_scatter_Mg.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随震源深度的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震源深度列表\n",
    "value_list = data[2][:, 1, 0]\n",
    "value_list[value_list < 1.0] = 1.0\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(value_min, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Depth (km)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='upper left')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_Sa_scatter_Depth.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随震中距的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震中距列表\n",
    "value_list = data[2][:, 2, 0]\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(value_min, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Depth (km)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='upper left')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, '%s_Sa_scatter_Rup.%s' % (dataset_type, fig_fmt)), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 特定台站结果可视化"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 指定台站，获取这些台站的数据，并计算残差"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取指定台站的data, label, pred数据\n",
    "stations = ['IBRH13', 'IWTH08']\n",
    "data, label, pred = [], [], []\n",
    "events = []\n",
    "for i in range(len(dataset.train_data)):\n",
    "    data.append([])\n",
    "\n",
    "count = 0\n",
    "for i, idx in enumerate(dataset.train_idx):\n",
    "    if dataset.train_events[idx][:6] in stations:\n",
    "        count += 1\n",
    "        events.append(dataset.train_events[idx])\n",
    "        label.append(dataset.train_label[i, ...])\n",
    "        pred.append(results['train_pred'][i, ...])\n",
    "        for j in range(len(dataset.train_data)):\n",
    "            data[j].append(dataset.train_data[j][i, ...])\n",
    "print('Train dataset contains %d events' % count)\n",
    "\n",
    "count = 0\n",
    "for i, idx in enumerate(dataset.valid_idx):\n",
    "    if dataset.valid_events[idx][:6] in stations:\n",
    "        count += 1\n",
    "        events.append(dataset.valid_events[idx])\n",
    "        label.append(dataset.valid_label[i, ...])\n",
    "        pred.append(results['valid_pred'][i, ...])\n",
    "        for j in range(len(dataset.valid_data)):\n",
    "            data[j].append(dataset.valid_data[j][i, ...])\n",
    "print('valid dataset contains %d events' % count)\n",
    "\n",
    "count = 0\n",
    "for i, idx in enumerate(dataset.test_idx):\n",
    "    if dataset.test_events[idx][:6] in stations:\n",
    "        count += 1\n",
    "        events.append(dataset.test_events[idx])\n",
    "        label.append(dataset.test_label[i, ...])\n",
    "        pred.append(results['test_pred'][i, ...])\n",
    "        for j in range(len(dataset.test_data)):\n",
    "            data[j].append(dataset.test_data[j][i, ...])\n",
    "print('test dataset contains %d events' % count)\n",
    "        \n",
    "label = np.array(label)\n",
    "label = label * dataset.label_std + dataset.label_mean\n",
    "pred = np.array(pred)\n",
    "for i in range(len(data)):\n",
    "    data[i] = np.array(data[i])\n",
    "    data[i] = data[i] * dataset.data_std[i] + dataset.data_mean[i]\n",
    "\n",
    "# 计算残差\n",
    "period_sec = [0.01, 0.1, 0.2, 0.5, 1.0, 2.0]\n",
    "resi_type = 'abs'\n",
    "\n",
    "if resi_type == 'abs':\n",
    "    residual = pred - label\n",
    "elif resi_type == 'rel':\n",
    "    residual = (pred - label) / label\n",
    "elif resi_type == 'log':\n",
    "    residual = np.log(pred) - np.log(label)\n",
    "\n",
    "data_scale = 'log'\n",
    "if data_scale == 'log':\n",
    "    label = np.log(10) * label\n",
    "    pred = np.log(10) * pred\n",
    "    data[0] = np.log(10) * data[0]\n",
    "    residual = np.log(10) * residual"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 不同周期下预测和实际反应谱值对比图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    r = np.corrcoef(label_points, pred_points)[0, 1]\n",
    "    xmin = np.min(label_points)\n",
    "    xmax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=label_points, y=pred_points, label='T=%.2fs' % T)\n",
    "    plt.plot([xmin, xmax], [xmin, xmax], 'k', linewidth=3)\n",
    "    plt.text(xmin ** 0.9 * xmax ** 0.1, xmin ** 0.1 * xmax ** 0.9, 'r=%.1f%%' % (100 * r), size=28)\n",
    "    plt.xlim([0.9 * xmin, 1.1 * xmax])\n",
    "    plt.ylim([0.9 * xmin, 1.1 * xmax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Target (gal)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Prediction (gal)')\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "    plt.legend(loc='lower right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_Sa_T_compare.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 绘制残差分布图"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随周期的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(period)\n",
    "num_bar = 10\n",
    "t_err, mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar), np.zeros(num_bar)\n",
    "for i in range(num_bar):\n",
    "    t_err[i] = period[int((i + 0.5) * N / num_bar)]\n",
    "    mean_err[i] = np.mean(residual[:, int(i * N / num_bar) : int((i + 1) * N / num_bar), 0])\n",
    "    std_err[i] = np.std(residual[:, int(i * N / num_bar) : int((i + 1) * N / num_bar), 0])\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "sns.scatterplot(x=np.tile(period, residual.shape[0]), y=residual[..., 0].ravel(), s=50, color='silver')\n",
    "plt.errorbar(t_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "plt.ylim([-0.6 * np.max(np.abs(residual[..., 0].ravel())), 0.6 * np.max(np.abs(residual[..., 0].ravel()))])\n",
    "plt.xscale('log')\n",
    "plt.xlabel('T (s)')\n",
    "plt.ylabel('Residual')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_residual_T.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随地表PGA的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载加速度时程信息数据\n",
    "event_msg = np.load('event_msg.npy', allow_pickle=True).item()\n",
    "# 获取数据集中每个地震动对应的PGA数据\n",
    "PGA_list = []\n",
    "for i in range(label.shape[0]):\n",
    "    PGA_list.append(981 * np.sqrt(event_msg[events[i]][13] * event_msg[events[i]][39]))\n",
    "value_list = np.array(PGA_list)\n",
    "\n",
    "# 将所有地震的PGA按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "num_bar = 10\n",
    "value_step = (np.log(value_max) - np.log(value_min)) / num_bar\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = np.exp((np.arange(num_bar) + 0.5) * value_step + np.log(value_min))\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if np.log(value) >= np.log(value_min) + i * value_step and np.log(value) < np.log(value_min) + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随PGA的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err[:], mean_err[:], std_err[:], fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('PGA (gal)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='lower right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_residual_PGA.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随震级的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震级列表\n",
    "value_list = data[2][:, 0, 0]\n",
    "\n",
    "# 将所有地震对应的震级按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.round(np.min(value_list)) - 0.5\n",
    "value_max = np.round(np.max(value_list)) + 0.5\n",
    "value_step = 1\n",
    "num_bar = int((value_max - value_min) / value_step)\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = np.arange(int(value_min + 0.5), int(value_max + 0.5))\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if value >= value_min + i * value_step and value < value_min + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随震级的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([value_min, value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Mg.')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.legend(loc='lower right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_residual_Mg.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随震源深度的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震源深度列表\n",
    "value_list = data[2][:, 1, 0]\n",
    "value_list[value_list < 1.0] = 1.0\n",
    "\n",
    "# 将所有地震对应的震源深度按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "num_bar = 10\n",
    "value_step = (np.log(value_max) - np.log(value_min)) / num_bar\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = np.exp(np.log(value_min) + (np.arange(num_bar) + 0.5) * value_step)\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if np.log(value) >= np.log(value_min) + i * value_step and np.log(value) < np.log(value_min) + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随震源深度的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Depth (km)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='lower right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_residual_Depth.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 残差随震中距的分布图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震中距列表\n",
    "value_list = data[2][:, 2, 0]\n",
    "\n",
    "# 将所有地震对应的震中距按从小到大的顺序划分为10组，以便后面求解误差棒\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "num_bar = 10\n",
    "value_step = (np.log(value_max) - np.log(value_min)) / num_bar\n",
    "sort_idx = np.argsort(value_list)\n",
    "value_list = value_list[sort_idx]\n",
    "N = len(value_list)\n",
    "idx_err = []\n",
    "value_err = np.exp(np.log(value_min) + (np.arange(num_bar) + 0.5) * value_step)\n",
    "for i in range(num_bar):\n",
    "    idx = []\n",
    "    for j, value in enumerate(value_list):\n",
    "        if np.log(value) >= np.log(value_min) + i * value_step and np.log(value) < np.log(value_min) + (i + 1) * value_step:\n",
    "            idx.append(j)\n",
    "    idx_err.append(idx)\n",
    "\n",
    "# 绘制残差随震中距的分布图和分组误差棒\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    residual_points = residual[:, idx, 0]\n",
    "    residual_points = residual_points[sort_idx]\n",
    "\n",
    "    mean_err, std_err = np.zeros(num_bar), np.zeros(num_bar)\n",
    "    for j in range(num_bar):\n",
    "        mean_err[j] = np.mean(residual_points[idx_err[j]])\n",
    "        std_err[j] = np.std(residual_points[idx_err[j]])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=residual_points, label='T=%.2fs' % T, s=50, color='silver')\n",
    "    plt.errorbar(value_err, mean_err, std_err, fmt='s-', ecolor='k', ms=10, mfc='r', mec='k', color='k', elinewidth=2.5, capsize=5, capthick=3)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([-0.6 * np.max(np.abs(residual_points)), 0.6 * np.max(np.abs(residual_points))])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('R (km)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Residual')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='upper left')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_residual_Rup.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 绘制label和prediction的散点对比图"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随周期的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ymin = 0\n",
    "ymax = np.max([np.max(label[..., 0].ravel()), np.max(pred[..., 0].ravel())])\n",
    "plt.figure(figsize=(8, 6))\n",
    "sns.scatterplot(x=np.tile(period, label.shape[0]), y=label[..., 0].ravel(), s=20, label='Target')\n",
    "sns.scatterplot(x=np.tile(period, label.shape[0]), y=pred[..., 0].ravel(), s=20, label='Prediction')\n",
    "# plt.ylim([ymin, 0.2 * ymax])\n",
    "plt.xscale('log')\n",
    "plt.xlabel('T (s)')\n",
    "plt.ylabel('Sa (gal)')\n",
    "plt.legend()\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_Sa_scatter_T.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随地表PGA的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载加速度时程信息数据\n",
    "event_msg = np.load('event_msg.npy', allow_pickle=True).item()\n",
    "# 获取数据集中每个地震动对应的PGA数据\n",
    "PGA_list = []\n",
    "for i in range(label.shape[0]):\n",
    "    if dataset_type == 'train':\n",
    "        name = dataset.train_events[dataset.train_idx[i]]\n",
    "    elif dataset_type == 'valid':\n",
    "        name = dataset.valid_events[dataset.valid_idx[i]]\n",
    "    elif dataset_type == 'test':\n",
    "        name = dataset.test_events[dataset.test_idx[i]]\n",
    "    PGA_list.append(981 * np.sqrt(event_msg[name][13] * event_msg[name][39]))\n",
    "value_list = np.array(PGA_list)\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(value_min ** 0.3 * value_max ** 0.7, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('PGA (gal)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='upper right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_Sa_scatter_PGA.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随震级的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震级列表\n",
    "value_list = data[2][:, 0, 0]\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(0.8 * value_max, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Mg.')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.legend(loc='upper right')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_Sa_scatter_Mg.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随震源深度的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震源深度列表\n",
    "value_list = data[2][:, 1, 0]\n",
    "value_list[value_list < 1.0] = 1.0\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(value_min, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Depth (km)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='upper left')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_Sa_scatter_Depth.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 随震中距的变化图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 获取震中距列表\n",
    "value_list = data[2][:, 2, 0]\n",
    "\n",
    "value_min = np.min(value_list)\n",
    "value_max = np.max(value_list)\n",
    "\n",
    "plt.figure(figsize=(18, 9))\n",
    "plt.subplots_adjust(left=0,bottom=0,top=1,right=1,hspace=0.2,wspace=0.2)\n",
    "for i, T in enumerate(period_sec):\n",
    "    idx = np.argmin(np.abs(period - T))\n",
    "    label_points = label[:, idx, 0]\n",
    "    pred_points = pred[:, idx, 0]\n",
    "    ymin = 0\n",
    "    ymax = np.max([np.max(label_points), np.max(pred_points)])\n",
    "\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    sns.scatterplot(x=value_list, y=label_points, label='Target', s=50)\n",
    "    sns.scatterplot(x=value_list, y=pred_points, label='Prediction', s=50)\n",
    "    plt.text(value_min, 0.12 * ymax, 'T=%.2fs' % T, size=24)\n",
    "    plt.xlim([0.8 * value_min, 1.1 * value_max])\n",
    "    plt.ylim([ymin, 0.2 * ymax])\n",
    "    if i >= 3:\n",
    "        plt.xlabel('Depth (km)')\n",
    "    if i in [0, 3]:\n",
    "        plt.ylabel('Sa (gal)')\n",
    "    plt.xscale('log')\n",
    "    plt.legend(loc='upper left')\n",
    "plt.savefig(os.path.join(model_path, result_path, figure_path, 'Stations_Sa_scatter_Rup.%s' % fig_fmt), dpi=300, bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}