FeatureSelection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "__author__ = 'Zhanwei Liu'\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from sklearn.preprocessing import MinMaxScaler,StandardScaler\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.model_selection import TimeSeriesSplit\n",
    "from sklearn.metrics import  mean_absolute_error\n",
    "from sklearn.model_selection import cross_validate\n",
    "from sklearn.metrics import make_scorer\n",
    "\n",
    "from statsmodels.tsa.stattools import acf\n",
    "from statsmodels.tsa.stattools import pacf\n",
    "from statsmodels.tsa.stattools import ccf\n",
    "\n",
    "from minepy import MINE\n",
    "\n",
    "# Custom module\n",
    "# from imp import reload\n",
    "# reload(MyPrepaData)\n",
    "from MyEvalCriteria import mae,rmse,corr,KGE,ia,BHV\n",
    "import MyPrepaData\n",
    "\n",
    "# Plot setting\n",
    "import seaborn as sns\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "rcParams['font.sans-serif'] = 'Times New Roman'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparation data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2011-2014 训练集；2015-2016 验证集；2017-2018 测试集。\n",
    "data = pd.read_csv('2011-2018再分析和流量数据.csv',index_col=0)\n",
    "data_train,data_test = data[:'20151231'],data['20151231':]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calc MIC\n",
    "mine = MINE(alpha=0.6, c=15,est='mic_approx')\n",
    "def mic(a,b):\n",
    "    mine.compute_score(a,b)\n",
    "    return mine.mic()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2011-01-01    420\n",
       "2011-01-02    415\n",
       "2011-01-03    435\n",
       "2011-01-04    415\n",
       "2011-01-05    433\n",
       "             ... \n",
       "2015-12-27    411\n",
       "2015-12-28    284\n",
       "2015-12-29    356\n",
       "2015-12-30    292\n",
       "2015-12-31    332\n",
       "Name: Q, Length: 1826, dtype: int64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_train.Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Determining the lag of runoff and rainfall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAChCAYAAAA2oJNvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAY+0lEQVR4nO3deXxV5Z3H8c8vYQlbQoAEwhKirILRChlZlKWKCy7MtI4LDHW0ReZldZQBx1q1LUyt1pla29dL6zR2lNqiLehQKigyouCG1qCyyhKEBAgQlhDWIMszf5xLSELg3Bvunu/79crrbs9zzu/k3uf+7nPOeZ5jzjlERETOJiXWAYiISPxTshAREV9KFiIi4kvJQkREfClZiIiILyULERHxpWQhIiK+wpYszKyJmT1tZj0Cj5ub2QQzG2dmvcK1HpFEYWYDzex6M7unxnM5ZjYnlnGJNESTcC3IOXfMzCoBCzw1HNgDLAIeBh6oWd7MJgITAVq1ajWwb9++4QpFJChLly7d5ZzLiuAqJjjn7jazyWZ2IbAN6AbsP1MFtQuJtTO1i7Ali3o0B44Cx4HDdV90zhUChQAFBQWuqKgogqGInM7MSiK8itTA7ck2MAbIAPqa2RDn3JK6FdQuJNbO1C7ClizMLBXoCOSY2eXAK8CdwGjghXCtRySBTDez0cAhYJhz7kUAMyuoL1GIxLNw7oY6DtwdePh+4Pa/z3W5eQ/NC6rcpp9ff66rEgkr59xHZ3h+fLRjETlXOhtKRER8KVmIiIgvJQsREfGlZCFJpbS0lGPHjgFQWVlJZWVljCMSib2XXnqJ559/HoDjx4+zdevWkJehZCFJY/369WzYsIEmTbzzNjIyMpg/fz5Hjx6NcWQisbN+/Xq6deuGmTcELjU1lY0bN1JSEtqZ40oWkjSee+45WrZsyezZs7n33nsByM3NZebMmTGOTCR2nnvuOZxzrFq1ittvv51PPvmEoUOH8otf/CKk5ShZSNLYtWsXubm5tG3btrqb3b59ezZt2hTbwERiaNeuXaSkpNC/f3+mTJnC4sWLSUlJCXkXbSRHcItEVZcuXXj//fcpLS3FzCgvL6esrIz8/PxYhyYSM126dCEzM5P169eTlpbG2LFjOXr0KDk5OSEtR8lCksaDDz7IsmXLuOWWW6qfe/PNNxk/XmPgpPE62S6efPLJ6ufeeust7rvvvpCWo2QhSSMzM5MBAwZw/PhxUlNT2b9/P2PHjiU1NdW/skiSqtsuTpw4waWXXkpmZmZIy1GykKSSnp5efb9NmzYxjEQkftRsFykpKSEnCtABbhERCYJ6FiIJSBNsSrSpZyEiIr6ULERExJeShYiI+NIxC5FGINRjHDomInWF87KqA4FOQJ5z7lkz6wtcAVQC851zu8O1LhERia5w9iwmOOfuNrPJZnYhUAF0Br4A9oZxPSISZ9QTSX7hPGZxcpjsceAw0Bz4M3A90K1uYTObaGZFZla0c+fOMIYhIiLhFs5kMd3MRgOHgGFAa2AQsKG+9TjnCp1zBc65gqysrDCGISIi4Ra23VDOuY/qeXp5uJYvkmjqOY6XAdwEXO2cuy220SU/7RoLL50NJRI5tY7jOedWmtlc4IwT85jZRGAieBduSmbBfJnX/CKP9Jd/MiSXSG6DkoVI5NQ9jgdwEOhgZr2dc+vqVnDOFQKFAAUFBS4qUUqDRCO5hJpQI0nJQiRyah3HM7NU4BvASmBzTCMTCZGShUiEnOE43mm9CWkcEn03l6b7EBERX0oWIiLiS8lCRER8KVmIiIgvJQsREfGlZCEiIr6ULERExJeShYiI+FKyEBERX0oWIiLiS8lCRER8KVmIiIgvJQsREfEVtlln614VLPDctcBOYJVzripc6xIRkegKZ89ignNuHtDczC40s6HAELzrcasHIyKSwML5JV73qmCDgFnAbmB03cJmNtHMisysaOfOnWEMQ0REwi2cyaLWVcGA2cBIoBXwXt3CzrlC51yBc64gKysrjGGIiEi4he2YxRmuCvZMuJYvIiKxo2MJIiLiS8lCRER8KVmIiIgvJQuRCDGzgWZ2vZndE3jc08x+aGZ/NbOOsY5PJBRKFiKRU2vsEVDinHsCWAA0i21oIqFRshCJnFpjj5xzR82sE/C5c25zfRU0/kji1VmThZm1jVYgIkmo1tgjMxsA/Ay4wswuq6+Cxh9JvPIbZzHWzPYAC5xzFdEISCRZnGHs0feiHohIGPjthloOlDvnKswsOxoBiYhI/PHrWUwF9pjZd4FueNN3iDRqgbObPgPWqMctjYVfsvg3YL1z7oiZpUcjIJEEsAzIcs4tMbMWzrnDsQ5IJNL8ksXlzrmVgft9gb9FOB6RRPAokGJmN+Ndw+WqGMcjEnF+yeLLGvczIhmISAKZ6JwrBTCz82IdjEg0+B3g7m9mKWaWAgyORkAiCWCcmaUF7mfGNBKRKPHrWbwBzAWaAj+IfDgiCeHDGpcJbh3TSESi5Kw9C+fcJufcdc65qwCdOivi6VLjfv+YRSESRb4XPzKzi4D7gTRgfsQjEol/vc1sAXAE+FOsgxGJhrMmCzP7FbAD+KNz7t3ohCQS91YB/4U395N6FtIo+B3g/gHwOXCe3zxRdadjDjyXY2ZzwhCnSDzZ7Zw77Jw7ALSLdTAi0eCXLG4FNgOvAN/3KVtrOmYza4836nt/fYU1u6YksIvMrHPgmhRDYx2MSDT4JYu2QO/ACNUnfMrWmo4ZGIPXkPqa2ZC6hTW7piSwGXizx/4SeCHGsYhEhd8B7lygjZllAr2AH56lbK3pmJ1zLwKYWYFzbklYohWJA8653cCdsY5DJJr8ksVyoCSYBZ1hOmacc+NDDUpEROLLWZOFc+6lGg8XRzgWERGJU7qsqoiI+FKyEBERX0oWIiJnkXriONet+YAfLXyeb274tPr5f/nkVQaVrmjwcnMrtvHTBb8Junza0Sp+M/vxBq/vXPlO9yEiEq+uW/MBB5q1IH97Mc8OvZUriv8GD70PY8Yw/vM3+DqlCR91v4jDTdPotbuUj3MvCnkdfcs3kl51gJ9eeVet5/emtcFwDY69NDOH1kcOBV2+qmkaR5o0a/D6zpWShYgkrGvWLWHyDZMZWrqcrAMV5G8vhnatoE8flnfqSauvq9jfvBWXlSyjokU6nfeVU5aeDc7xrVXv0vT4Mfi8M6xcyYAtO7i85AsWnT+Q25YtoKpJM+b1vZyeuzfTfe92ulSWc8+SmXD3hfDuu/TaVUpJZg5Xr1tC1sEK1nfI5Zp1S1jRqSfXrFvCpBsfYPTaD+GpNXDbbYxZvYi2h/fzabf+fJl9PgDtD1Uy+b0/sLtVWw41TWNLRkd4fBlcdRX87nf8eOl25vW9nJz9u9jWpgNZB2N3FV/thhKJkDNMgXOrmd0Rw7CSygsFY7jhy/fIqyhjd8t0fn35OHjwQZgxg+U5vVnS/SIuK1nGcUthe5v2jFn9HgDtDu+jbdV+Zl58NfTrBwsX8lnXC+hauYM9LTPY3TKDBb0H02P3Fkoyc9iUmcPWjGxaHK2CWbNg4EDWd8gFYORXRcy45DrWZOWxv3lLPu/ch4PNWnBB+UYONmsBU6ZAkyb02rWZlwbeyNb0UxN4726ZwS+HjeeSsjWszj6P/c1bwurVkJUFWVnVMXxzw6cs7dqPna1id/kUJQuRyKk1BU7guR2oRx82O1q353DTNP5r+O2kHfuaO4r+CsuXw803A3Bl8Se806OANdl55G8vZnNGRwD2tEin165Svr9kJlRWQs+eDP9qKR92v5jsA3tof2gvmYf20eHQXjrv20nnfTvpcLCCtlX7YdgweOYZ8iq20fbwfo6lNOHHbxfSZV852QcqyD6wh4yqA2xr0z7Qs3gKWrbk/D1b+PfFvye96kB1/IZj2KbPmZV/Fdet/ZAuleVw6BCUlEB5eXUMWzI6ctsX82lz5CBtjhyMyf9aH1qRyKk7BY4vM5sITATIzc2NUFjJY3t6B7and6h+PL1gDFO/fX3144U9BwFQktmZkszOpyqa8fC1/wrAg9nZ8OijvPfQvOqXP+tyQb3ru/PmaWwaNAgGDeLJQPn5fS6rfv3ha+8F4NNu3m+DyTdM4dtTvHju+YfTJ8CYdOO/V9//MO8b3p2f/8y7HTGCN2vEBPCnb1xbb1zRoJ6FSOTUmgIn8FxnIOdMFTRnmsQr9SxEIqS+KXCccy/HIhaRc6VkISKNQnrVAa5d+xHnVZQxp98I1mSfF+uQEoqShYgkrbSjVVxZ/CkX7tjAvrRWzO891DsDSkKmZCGnHDsGs2fDRx/BqFFwfeBA4X/+JwwaBCNGNGy5GzZ4Z4T8JsjRqocOwe23w6uvNmx90rgdPQoLF/Lg4ukcSW3G2z0vZV7fy8Es1pElNCWLRDFrFqSnQ1ERPPIIzJ0LH3wAY8bAsmXQrBlccQW0auWdpz1yZOjrWL4c9u6Fp5+u/Xy7duAaPlKVHj1g377gy7dsCWlpDV+fND7OMWDrGkZsXArTlsCoUfxi2Hc4kZLqX1eComSRKGbPhpdegnfege3bYelS2LMH+vSBpk3hwAHIyICFC6FDBygthdxc70v+j3+Er7+GAQNg5Urvy/vtt2H0aPjd76BFC++89C+/9HoBJSXw+OPw8MPw7rte8unZE/7yF2/d/fp58Qwc6N3OmOH1AnbsgNtug8WLvdiGDYOLL/bi37kTfvQjyM6G1q0hLw+WLKkeqVodw5Yt0LWrtx4RH312buKq9R/T8mgVn3W+gGeH3Mr9j/0DACfmz/OpLaFQskgUkybBn/8M69d7ozt/8hMoLva+qO+7zyszaxakpkKXLvDKK/CDH8CuXd4X9/33w5Ej8Otfw3e+431BB0aJMmoUrFnjJQSA7t3h4EFveddc4+2eAnjjDSgshIoKL2kNHuwlnWXLoE0bGD/eSxirV8Njj3nlTsrKgv/4D6/MlCleElu9GsaNqx3D4sVeUuzUKbr/X0kYXfdu57q1H5J5eD9rsrrzQsHfc6hZi1iHlfTCkizMbCDQCchzzj1rZj2Bm4EhwF3OuR3hWE+j1rkzlJV5v/gPH4YXXvB+gQdGqvL6694xhrIy+PRTOC9wpkeHDrBqlVdvwgQvIbz1Flx5JWzbBuXlXkLZscProZSWevdP9gyeeQYyM72/pk29pHXnnV7dbdu8hNC1q1euuBgmToS1a71eyV13efXASw7/939e3VdfhYKCWiNVq2PIy4Pnn/dG1VZWer0laZQKtqyC4j7eZ7a8nO98NpfO+3axJSObWfmjqGipz0Y0hatnMcE5d7eZTQ5Ma7DWOfeEmd0L1DtNokaqhqhrV+/vpJO9iZNuvNG77dnzVA8BvIN6hYWnHj/6aO16Q4bUv7433vBuBw069dxNN526/9vferfDAmPN/vCHU6/NmnX68mbMOHV/1Cjv9tvf9m7rO3B+112nPyeNxiVb1zB84+fwQbbXC87OZmHPS71JACUmwjWCu9a0Bs65o2bWCfjcObe5vgoaqSoi9RmwZTW/mfME2Qf2eLtUH38cJk9WooixcPUsak1rYGYZwD3AJjPDOfdhmNYjIknui5ze3Dz253Tbt5PbNm/2dnW2bx/rsBq9sCSL+qY1AL4XjmVLmKxc6R2HGD481pGInNWJ1CZsycxhS2YOfPd6/woSFTobKpmdOAHz5sF770H//t5prZK0BpcuZ3DpCt7scxmDS1eQeXgflF0ChYVcu/ZrDjVNY0DZGv56wQiu2PA3mh0/xl/7jeAfV7zNik6B41xTP4U77uD7S2bydWoT3ulxKWO+XMxnnfvS8mgVfXdu4uWLr2XcsvlUtEiHFbnw2msMLm1O9oE9nL9na/Xr5a3bsbxTL0YVf8J7eQM4v2Krt/yJE6GwkH9cUclXmV0Yvukz3u45iIu2ryf7wJ7q+l+16wKLWsGiRXDTTfzz0tfJPLyv+vU1WXn1bhOb+sP06VxZfAKA/O3FvJo/ijGrF9e7TSdjmvTBDCpapPNxbj6j137Ix7n59W4TSzvB668zYEtrzq/YStfKHdWvb8noeNo21Vz+V+26UN663WnvU81t4q0m3mnlY8cy4W//W+/7VHObeGoN3HADvPIKw79KOe19+jg3v0FXCKzL3LkMtgqTgoICV1RUVO9reQ8Fd670pp/rF0i1ffu80083b4brrvN6Exq9ehozW+qcK4h1HGcSznYRb+WDrRNq+XOJKd7KB1vnXP5H9TlTu1DPIpkUF3sD8FJSvOky8vJiHZGIJAkli0TnnDdq+623vJHZDzzgjZAWEQkjJYtEdeiQN3Zh7VpvgN2TT3o9ChGRCFCySDSbN8Pvf+9N3TFunAaviUhUKFkkAue8acPnzPGm/bjnnlPTaIiIRIGSRTz7+mtv8sAvvoDLLvNGsjbRWyYi0advnnhx5Ag0b+7dLy+HF1/0ri1xyy3eLLEiIjGkZBEPCgshP9+74M/Mmd4upjvugI4dYx2ZiAigZBF7b77p7WZauxYuuQSmTfOueicJr56p+zOBvweaAq855/bENECREChZNEBYR5U/9JB3xbkOHbwr2ylRJJO6U/cPAZYGXhsHPBO70ERCo+k+GiAeY5LQRXq6DzMrdM5NNLP7gbnAaGAx3qUBCpxz/1NPnerrvGRkZAycNGlSpMITqde0adOSZ7qPs02YRn6+N4I5MBEXc+eeGpMwfbr3Kx68a1jfcQe8/DKPvb3Bd8K0ab+8B157DUaOZMzqRb4TpnWt3HEqprw86N0bFizwLlJUVORdZS4woRq9e3unxAYmTGPRIti9+9TrDdgmmjevnlyMIUO8a3SvWHFqme3bw8iR1dtEWRmsW3fq9Zwc72p2r78OV1/tvbZp06nX43GbRo70/uJHran7gT8BgatUMbu+Cs65QqAQvB9RU6dOjUKYIqdMmzat3ufVswhyHeGerCvRJOP/KJEnEhSJFE0kmMSiMZuliDRuShZRkAi/skVEziZsyUKnCYqIJK+wHbMws+dOniYILKD2aYJDnXPP1Cmvsz6k2q/eXhdUuUmjejeofH3OdNZHvNAxC4mFaByzSA3cHgcOA82Bo3inCR6uW1hnfUhNob7906uC27U3deqZd+2d6awPETldOJNFyKcJiohIYghbsnDOfVTP0y+Ga/kiNelkAJHo0qXVRETEV2KeOrtokUY7a5sSYQS3SNKI+xHcIpGSdCO49SNK2xSGH1FnahdKFtJoJV2yEAmDM7ULHbMQERFfShYiIuJLyUJERHzFxTELM9sJlIRQpQOwK4Llo7GORC8fjXVEunx351xWCOWjKg7bhT5TiVe+IXXqbxfOuYT7A4oiWT4a60j08vEYU0O2IZn+4u3/2xg/U/FWvqF16vvTbigREfGlZCEiIr4SNVkURrh8NNaR6OWjsY5obEMyibf/b2P8TMVb+YbWOU1cHOAWEZH4lqg9CxERiSIlCxER8ZVwycLMmpjZ02bWI8jyI8zsZ2b2hyDL9zSzSWb2WAgx5ZvZ/wRZtkUgnh+ZmQVZJ9XMbjGz/kGW/5aZPWFmc8ysaxDlx5rZVWb2vSCXf7OZXWNmjwRRtvr9MrPxZvZNM7shyPId/N63OuXvNbMXzeyfgtmOZBJKuwi1TQTqxFW7iHSbCNRRu6gh4ZKFc+4YUAkE9UULfOScewQoC3L5xcA7wJpgyptZb+Ag3mVkg3ER0BEodsEfMHoA77K1Q4Mp7Jyb7Zz7ITDHObcliCq7ge7AsiDj+RawEEg3s04+sZx8v1oCA5xz7wI3BVHenHO78Hmfa5RPAf4bmAJ0CXI7kkaI7SKkNhFYfry1i0i3CVC7qCXhkkWonHNHzexC4IUQqm0Hrg6y7K3AGKCvmXULIp5PnHMTCPJDHtDbOfcKMNjM2gdTIfBhDfrLAJgL3GlmwXwmfoW33b2AnUEuPwXv+uwAB0KIK1gusPyhwNMRWH7SaGCbgPhqF9FoE6B2US3hLn5kZql4v0BygOIgyo8HLgX2mtmzzrkdPuVHAc2AOcHE45z7aaBegXNucxDx3ACk4X0IgzXbzG4EVjjndgdZ5yrg5SDLZgP9gE1Bli8DOgMPO+eOn61gjfcrHVhtZsOBeUGUzzGz3UC2mbV2ztXbkOp8Hv4N2IH3a/DZILclKYTSLkJtE4E68dYuIt0mQO2i9jJ16qyIiPhJ+t1QIiJy7pQsRETEl5KFiIj4UrIQERFfShaNgJk9a2Y9Yx2HSDxRuwiNkkWMmNmdZvZTMwvq9OXAiMwHzOzJwOMeZvbdIFcX7DnfIjGldhG/lCxiZyNQEhhpiZkNCkxJ8Ezg8VQz+xczexqqR2QWASsDUy6UABvN7Ldm1sfM3jCz9MB0BhPM7Ckze8TMTo4KvdrM/mJmGWZ2u5ndb2b9zGyemd1tZmnR/xeInEbtIk4pWcSBwACaUmAvp4bktwZWAk/VKf4K8CXwE7zRmWV4IzTLgUPAfmA60B74M9A3UG8B8G7gcTPn3K+BbUCFc+4551xVJLZNpKHULuKLkkXs9AIGmNlo4F+BYcBAwJlZNtAUuAL4uxp1egDnOedmAOvwGsNq4HagFV6D6gR0A9oF7nfE+7U1HKjCa1A3mdmP8RpeZzPLjOymigRN7SJOaQR3HDKztsAEYAtQ5Zz7S4xDEok5tYvYUrIQERFf2g0lIiK+lCxERMSXkoWIiPhSshAREV9KFiIi4kvJQkREfP0/vCNR5MXryUwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 396.72x162 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def mypacf(x, nlags=40, method='ols-unbiased', alpha=None):\n",
    "    Qpacf = pacf(x=x, nlags=nlags, method=method, alpha=alpha)\n",
    "    if alpha is not None:\n",
    "        '''\n",
    "        Calculation method of confidence interval：\n",
    "        varacf = 1. / len(x)  # the length of series\n",
    "        # ppf: Percent point function (inverse of cumulative distribution function(cdf))\n",
    "        interval = stats.norm.ppf(1. - alpha / 2.) * np.sqrt(varacf) (while alpha=0.05, ≈ 1.96/√N)\n",
    "        \n",
    "        ref1: https://www.itl.nist.gov/div898/handbook/pmc/section4/pmc4463.htm\n",
    "        ref2: https://www.statsmodels.org/stable/_modules/statsmodels/tsa/stattools.html#pacf\n",
    "        '''\n",
    "        interval = abs((Qpacf[1][:,0]- Qpacf[0])[1])\n",
    "        return Qpacf[0],[interval for i in range(nlags+1)],[-interval for i in range(nlags+1)]\n",
    "    else:\n",
    "        return Qpacf[0]\n",
    "def myccf(x, y, nlags=50,unbiased=True, alpha=None):\n",
    "    cc = ccf(x=x,y=y,unbiased=True)\n",
    "    if alpha is not None:\n",
    "        varccf = 1. / len(x)\n",
    "        interval = norm.ppf(1. - alpha / 2.) * np.sqrt(varccf)\n",
    "        return cc[:nlags+1],[-varccf + interval for i in range(nlags+1)],[-varccf - interval for i in range(nlags+1)]\n",
    "    else:\n",
    "        return cc[:nlags]\n",
    "\n",
    "def myplot_pacf(x, nlags=40, method='ywunbiased', alpha=None):\n",
    "    if alpha is not None:\n",
    "        Qpacf,up,bottom = mypacf(x=x, nlags=nlags, method=method, alpha=alpha)\n",
    "        plt.plot(Qpacf)\n",
    "        plt.plot(up)\n",
    "        plt.plot(bottom)\n",
    "    else:\n",
    "        Qpacf = mypacf(x=x, nlags=nlags, method=method, alpha=alpha)\n",
    "        plt.plot(Qpacf)\n",
    "\n",
    "def myplot_ccf(x, y, nlags=50,unbiased=True, alpha=None):\n",
    "    if alpha is not None:\n",
    "        Qccf,up,bottom = myccf(x=x,y=y,nlags=nlags,unbiased=unbiased,alpha=alpha)\n",
    "        plt.plot(Qccf)\n",
    "        plt.plot(up)\n",
    "        plt.plot(bottom)\n",
    "    else:\n",
    "        Qccf = myccf(x=x,y=y,nlags=nlags,unbiased=unbiased,alpha=alpha)\n",
    "        plt.plot(Qccf)\n",
    "Qpacf,up,bottom = mypacf(x=data_train.Q.values,nlags=12,alpha=0.05)\n",
    "Qccf,upccf,bottomccf = myccf(x=data_train.real_p.values,y=data_train.Q.values,alpha=0.05,nlags=12)\n",
    "'''\n",
    "Plot Partial autocorrelation function（PACF）and Cross correlation function（CCF）\n",
    "'''\n",
    "fig = plt.figure(figsize=(5.51,2.25))\n",
    "rcParams.update({'font.size': 7})\n",
    "plt.subplot(121)\n",
    "plt.bar(range(1,len(list(Qpacf))),list(Qpacf)[1:])\n",
    "plt.plot(range(1,len(up)),up[1:],ls='--',c='r',lw=0.5)\n",
    "plt.plot(range(1,len(up)),bottom[1:],ls='--',c='r',lw=0.5)\n",
    "plt.hlines(y=0,xmin=-1,xmax=100,lw=0.5)\n",
    "plt.xlim(0,13)\n",
    "plt.xticks(range(1,13,1))\n",
    "plt.xlabel('Lag Number')\n",
    "plt.ylabel('PACF')\n",
    "plt.ylim(top=1)\n",
    "plt.yticks(np.linspace(0,1,6))\n",
    "plt.annotate(\"95% confidence band\",xy=(5,0.05),xytext=(6,0.2),color=\"r\",\n",
    "arrowprops=dict(arrowstyle=\"->\",connectionstyle=\"arc3\",color=\"r\",linewidth=0.5))\n",
    "plt.text(s='(a)',x=12,y=0.9)\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.bar(range(1,len(list(Qccf))),list(Qccf)[1:])\n",
    "plt.plot(range(1,len(upccf)),upccf[1:],ls='--',c='r',lw=0.5)\n",
    "plt.plot(range(1,len(upccf)),bottomccf[1:],ls='--',c='r',lw=0.5)\n",
    "plt.hlines(y=0,xmin=-1,xmax=100,lw=0.5)\n",
    "plt.xlim(0,13)\n",
    "plt.xticks(range(1,13,1))\n",
    "plt.yticks(np.linspace(0,0.4,5),['%.1f'%i for i in np.linspace(0,0.4,5)])\n",
    "plt.annotate(\"95% confidence band\",xy=(5,0.05),xytext=(6,0.1),color=\"r\",\n",
    "arrowprops=dict(arrowstyle=\"->\",connectionstyle=\"arc3\",color=\"r\",linewidth=0.5))\n",
    "plt.xlabel('Lag Number')\n",
    "plt.ylabel('CCF')\n",
    "fig.tight_layout()\n",
    "plt.subplots_adjust(wspace=0.25)\n",
    "plt.text(s='(b)',x=12,y=0.36)\n",
    "plt.savefig('PACF-CCF.tif',dpi=300,bbox_inches='tight',pad_inches=0.05)\n",
    "print(end='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\software\\Anaconda\\lib\\site-packages\\statsmodels\\tsa\\stattools.py:572: FutureWarning: fft=True will become the default in a future version of statsmodels. To suppress this warning, explicitly set fft=False.\n",
      "  FutureWarning\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALQAAACiCAYAAADoQue0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAPZUlEQVR4nO3deXSUVZrH8e8TIAkJJISQQAKGCAGiGBSSNiwNLqMoCtjdLt06HJeRQx/bXtB2PC4zrU6PC2N309NHRqVn1BkPMh7boRWZcRjcQUQTWlYDsiQxhCUkIUBI2Hzmj7cIJCRVFfJWlpvnc05OUvXeuve+lV/dc6veeu8rqooxrojq6A4Y4ycLtHGKBdo4xQJtnGKBNk6xQBunWKCNU3wLtIj0FJH5IjI8cDtGRGaLyG0iMsKvdowJpqdfFanqCRGpASRw1xSgCvgQeAR44MzyIjIHmAMQHx+fm52d7VdXjOMKCwv3q2pKc9t8C3QzYoDjwEmgrulGVV0ILATIy8vTgoKCCHbFuERESlra5ueUowcwEEgTkTuB/wMGA9OAl/xqx5hg/JxynATuCdz8JPD7hbbWm/nQsqDbi5+5vq1NGIfYpxzGKRZo45RIvilsV6GmJmDTk+7ARmjjFGdG6Naw0dxdNkIbp3TLEbo17GPDrsVGaOMUG6F9ZKN5x7NAdwB7Uxo5NuUwTrERupOz0bx1LNAOsfBboLstV9/AWqBNSF0p/BZo45vOMOWxQJsOEanw28d2xim+jdAikgsMAjJVdYGIZANXAjXAu6pa6VdbxrTEzxF6tqouA2JE5CLgEJAO7AQONC0sInNEpEBECioqKnzshunO/Ax0j8DvU8sWxACvA9cD5zUtrKoLVTVPVfNSUppdYsGYVvMz0K+IyDTgCDAZ6APkA9t9bseYFvm5jMGnzdy93q/6jQmHjZzGKRZo4xQLtHGKBdo4xQJtnGKBNk6xQBunWKCNUyzQxikWaOMUC7RxigXaOMUCbZxigTZOsUAbp1igjVMs0MYpETvrO3DftUAFsElV6/1qy5iWROysbxGZCEzAO7/wrHbsrG8TCZE86zsfeAOoxLvedyN21reJhEie9b0EuByIBz72sR1jWhTps76f86t+Y8Jhn3IYp1igjVMs0MYpFmjjFAu0cYoF2jjFAm2cYoE2TgkaaBEZ3V4dMcYPoUbo0SJyl4gktUtvjGmjUIHeDRQDNSIyMvLdMaZtQn2X42FgD3A7MBC4LuI9MqYNQgX6aaBAVetEZFB7dMiYtgg15UhR1brA3xmR7owxbRUq0FVn/J0cyY4Y44dQgR4jIukiMhDvdCpjOrVQc+hFwG+AaGBl5LtjTNsEHaED1+degndV2LOuBnsmEckVketF5N4z7ksTkbd86akxYQh1pPB5vKUJfqOqj4Soq+lZ38l4L4JD/nTVmNBCjdD3AJvxjhhmhair6VnfM4GJQLaInDX/tmUMTCSE/HKSqq5U1X8Drg5RtNFZ36r6sqr+HihS1dXN1GvLGBjfhX3Wt6o+H2J7c2d9o6qzWtspY86VfX3UOMUCbZxigTZOsUAbp1igjVMs0MYpFmjjFAu0cYoF2jjFAm2cYoE2TrFAG6dYoI1TLNDGKRbobqTHtyfhjTfgvvu4YvsXDff/eM2fyC/dcO4Vb98OP/lJ+OWPHOFfljx17u0F4dtVsEzbXFe0ksPRvcnZs40FE38I77wDK1cyrmwAF1bs5FhUTz4dOoa6XrGMqCzls4wxrW4je99OOBAH8+fzwUPLGu4/ENsXQc+988OHw8GDkBBm+bg4jvaMPvf2grBAdxLXbF3N/dPvZ2LpelIOV0PhdqiqYkfyxZzo0YP4Y/UciolnUsk6qnsnkH5wH+UJqaAKr74Kx44xem8tIytKKOmXxndLvuTDYbn8aN1y6ntGsyz7u2RVfgPb46CkhKfefY4FE25hQul6RuwvpSQpjalbV5NSW83XAzLgvhWQmwtLlhBz/iymbVnFgNpqll4whfxvNtKv7hBfnDear1KHeTtQUcH9O16lMr4fR3rFUpY4kLHlRazMvKRRH3j9MAwZQkptdUSeR5tydBIv5c1k+lcfk1ldTmVcAjz2GDz4IN/b9CHr00ayeugYJpWs46REsadvMjM3e9cy7V93EKqq4O67+To5g0kl61g75AKG1OylKi6RyrhElo8cz/DKMkqS0iArC4YOpffxeq4rWsmGQVlegIHLdxSwaOx1FKVkQmIijB8Pfftywb6d1Eb35l8v/QEno3owYv83/EfuDHYlpJ7egZQUfjd5FmPLi9icej6HYuIYsb/0rD6wbBlMmkRFfGQWtPUl0E2XMBCRLBF5WETeDixSY0LY2yeZul6xPDvldmJPHIM//AHWr/dGNeCvtq3h/eF5FKVmkrNnG98kek9rVe8E2LQJnnqKvkdrKU5KY8qOQlYNvZjUw1UkHzlA0pGDDDhygPSDFVBaCnv30q/eG2HvKHyHzOrd9Ks7xImonvxqxUIGH9wHu3d7P9XV7O6bzLQtq5j9+X9R1zOGYVVl/O1H/05C/eHTO6DK5OK/8EbO1Vy3ZRWDa/bR+8RRBtfsa9QHMjPhj3+k79Fa+h6t9f159GvKMVtV7xGR+0XkImCLqj4tIj/FW6TmLCIyB5gDkJFhy+btSRjAnoQBp+/4+c8BqPjcm+u+l5UPQElSOiVJ6afLicDChQBUPrSM5yb+qFG9awdf0Oj27/7hegDuuvkJAL5MH9Ww7d1Rk04XfOZn3u/Jk9n70DLun/7Lhk33fu/hs3dg0SI+CczLV2VeAsD/jpoIwJqMnIZiDwban709nUjwa8rRaAkDVT0eWK30L6r6TXMPsLO+TST4NUI3WsJARBKBe4FiEUFVV/nUjjFB+RLoFpYwuNuPuk2Y6usZs3sr69O694UW7GO7Lq5f3UF49lk4cIDa6KDLD3YLFuguKqN6NzdufM87QPHwPBg0iO1nHCzprizQXczYXUVcs/VTyhIH8kL+TdRFx/LgILtayCkW6C4g6tuTXLXtc75Ttom16dn802V38G1Uj9AP7IYs0J1ZXR23fvkuQ6vLWTEinyevuNv73Nm0yALdCfU/UgPz5sHhw6w57yIWX3JtR3epy7BAdyLnV+3iBxvfpza6N/zdPEhNZYe90WsVC3RHUyVv12au/noNxUnpLJhwM/W9YnkoNTX0Y81ZLNAdJOrbk1yzdTXjyov4Ysho5tkbPV9YoNtJXtkmeL835Ocza+0yhtTsZfmICTx55eyO7ppTLNDtIO1gBTM2fwyrj8MHH/DJ+WMbf2PO+Ma+4N8OFi9+hIyaPRAXBzfcYGGOIBuh28E1dy8g9vhR1t0yGeLjgb0d3SVnWaDbwdGe0d53LgYP7uiuOK/LBXp86XrGl27gf0ZNYnzpBnj8C5gzh7krF1GUksmRXrGMKy/i7Qsu48rtnxN98gRvX3gZN21YAUu/BWDuysX8KecqZm7+iGM9evL+8EuZ+dVHrE3PJu54PdkVxVA+FhYu5I7CPXyWkcO0Lav4LCOH1MNVDKvaxWsXX8tt696FF8sgLw+WLmVcWR+GVe9iSM3ehu1liQPZkTSYKcVroXAQFBQwd+VHDdt39B/Mvj79G+1TUt3Bhvav3XIs6D5tGJTl7VdhIUNqMoLu02sXXwuPPw7JyYyqiGpxn/b16e/1delSmDqVmzasaHGfVmTlw4svwu7dpB4aFnSfGtrPyWHKjk2N9mnjoKxzOpO9KVFtw+nrPsnLy9OCgoJmt2WGOLBQ/Mz1YZWLVNlT5VpTtqP62pqynWm/mhKRQlXNa26bvSk0TvFtyiEiuXjXBc9U1QUikgTcAPQC3lTVKr/aMqYlvk05ROT5U2d+A8uBCUBhYPNEVX2uSfmGs74TExNz586d60s/jPueeOKJFqccfgZ6oarOEZFfAO8A04CP8KY1eYHrhTcr2BzamKaCzaH9/JSj0ZnfwH8CMwLblvjYjjEt8i3QLZz5/bJf9RsTDvuUwzilyx1Y4cMPvZ8bb/R+V1bCnDneclg5OdCnD6xeDbfe6i1Je/Qo3HYbvPKKt5omQGEh3HknvPYaxMTA9OmweDFMmACHD8OGDafrTE6Gyy+HN9/0fpeXw9atp7enpTUcWGHqVG9bcfHp7ZmZMHIkLF8OM2ZAQYG3Ztyp7SNHQnq67VNurteXNur0B1aMacoOrJhuwwJtnGKBNk6xQBundIo3hSJSAZSEWXwAsN/nspGoM1JlO7r9SJVtTZ1DVbX5RcVVtUv9AAV+l41EndbXyO1XsB+bchinWKCNU7pioBdGoGwk6oxU2Y5uP1JlW1NnizrFm0Jj/NIVR2hjWmSBNk7pUoEWkZ4iMl9Ehocod5mIPCkir4ZRZ5aIzBWRfwyjbI6ItHjmTZOyvQN9+HuRllcpF5EeInKLiIwOo87vi8jTIvKWiAwJUfZWEblaREJejUxEbhaRa0Tk0SBlGp57EZklIleIyPQwyg4I9n9oUvanIvKyiPx1qD63pEsFWlVPADVAqGXsP1XVR4HyMOrcBrwPFAUrJyIjgVogJrzeMgYYCGzT4G9UHsC7cOnEMPq6RFUfBt5S1bIQxSuBocC6MPr6feA9ICFwwdTm2j713McB41T1A+DGEGVFVfcT5P91Rtko4AXgl8A5r8jTpQIdLvWuZHsR8FKYD9kDTA1R5ofATCBbREJeP01V16jqbEIHdaSqLgbGi0hyqHoDgQv5Qg14B7hLREL9n3+Pt38jgIoQZaPwrhgMcDhYwVbSQL0TgfnnWkmX+oK/iPTAG/XSgG1Bys0CLgUOiMgCVW1xMTkRuQrveuRvBWtbVX8dKJ+nLVzuuUm904FYvFAFs0REZgAbVLUyVL3A1cBrYZRLBS4EisMoWw6kA4+o6snmCpzx3CcAm0VkCtDsajFn/p9EpBJIFZE+qnrWC6DJ//Q+vIX/hgILwuj32W3bx3bGJU5OOUz3ZYE2TrFAG6dYoI1TLNBdnIgsEJGsju5HZ2GBjgARuUtEfi0iYX0sGjha9oCIzAvcHi4ifxNmc6E+N+5WLNCRsRMoCRwFQ0TyA4etnwvcflxEfiwi86HhaFkBsDFwCL4E2CkiL4rIKBH5bxFJCBzyni0ivxWRR0Xk1JG6qSLyZxFJFJHbReQXInKhiCwTkXtEJLb9n4KOYYGOsMCBg1LgAKcP6fYBNgK/bVJ8MfAV8BjekbNyvKNn+/AWwTwEvAIkA68D2YHHLQc+CNyOVtV/BnYD1ar6vKrWR2LfOiMLdGSMAMYFVmP9Gd5qrLmAikgq3iLwVwLfOeMxw4HzVXURsBUvyJuB24F4vBfDIOA8oH/g74F4o/kUoB7vxXCjiPwK70WTHlh4vtuwI4XtTET6AbOBMqBeVf/cwV1yigXaOMWmHMYpFmjjFAu0cYoF2jjFAm2cYoE2Tvl/kgp+SVMTc8IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 180.72x162 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def myacf(x, nlags=40, method='ywunbiased', alpha=None):\n",
    "    Qacf = acf(x=x, nlags=nlags, alpha=alpha)\n",
    "    if alpha is not None:\n",
    "        '''\n",
    "        Calculation method of confidence interval：\n",
    "        varacf = 1. / len(x)  # the length of series\n",
    "        # ppf: Percent point function (inverse of cumulative distribution function(cdf))\n",
    "        interval = stats.norm.ppf(1. - alpha / 2.) * np.sqrt(varacf) (while alpha=0.05, ≈ 1.96/√N)\n",
    "        \n",
    "        ref1: https://www.itl.nist.gov/div898/handbook/pmc/section4/pmc4463.htm\n",
    "        ref2: https://www.statsmodels.org/stable/_modules/statsmodels/tsa/stattools.html#pacf\n",
    "        '''\n",
    "        interval = abs((Qacf[1][:,0]- Qacf[0])[1])\n",
    "        return Qacf[0],[interval for i in range(nlags+1)],[-interval for i in range(nlags+1)]\n",
    "    else:\n",
    "        return Qacf[0]\n",
    "Qacf,up,bottom = myacf(x=data_train.Q.values,nlags=12,alpha=0.05)\n",
    "fig = plt.figure(figsize=(2.51,2.25))\n",
    "rcParams.update({'font.size': 7})\n",
    "plt.bar(range(1,len(list(Qacf))),list(Qacf)[1:])\n",
    "plt.plot(range(1,len(up)),up[1:],ls='--',c='r',lw=0.5)\n",
    "plt.plot(range(1,len(up)),bottom[1:],ls='--',c='r',lw=0.5)\n",
    "plt.hlines(y=0,xmin=-1,xmax=100,lw=0.5)\n",
    "plt.xlim(0,13)\n",
    "plt.xticks(range(1,13,1))\n",
    "plt.xlabel('Lag Number')\n",
    "plt.ylabel('ACF')\n",
    "plt.ylim(top=1)\n",
    "plt.yticks(np.linspace(0,1,6))\n",
    "plt.annotate(\"95% confidence band\",xy=(5,0.05),xytext=(6,0.2),color=\"r\",\n",
    "arrowprops=dict(arrowstyle=\"->\",connectionstyle=\"arc3\",color=\"r\",linewidth=0.5))\n",
    "plt.savefig('ACF.tif',dpi=300,bbox_inches='tight',pad_inches=0.05)\n",
    "# plt.annotate(\"95% confidence band\",xy=(5,0.05),xytext=(6,0.2),color=\"r\",\n",
    "# arrowprops=dict(arrowstyle=\"->\",connectionstyle=\"arc3\",color=\"r\",linewidth=0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=4 face=\"宋体\">使用GBRT模型选择径流的滞后项与降雨的滞后项，首选选择一组默认参数，然后对各组特征进行计算选择误差较小的一组特征作为输入特征。2011-2014年为训练集，2015-2016年为验证集。根据PACF和CCF的分析，选择径流滞后1，4阶必选，降雨滞后1，2，3...20阶段备选,共21组备选输入结构:</font>\n",
    "- 1.$Q_{t-1}$, $Q_{t-4}$\n",
    "- 2.$Q_{t-1}$, $Q_{t-4}$, $R_{t-1}$\n",
    "- 3.$Q_{t-1}$, $Q_{t-4}$, $R_{t-1}$, $R_{t-2}$  \n",
    "- ...  \n",
    "- 12.$Q_{t-1}$, $Q_{t-4}$, $R_{t-1}$, $R_{t-2}$, ..., $R_{t-11}$\n",
    "- 13.$Q_{t-1}$, $Q_{t-4}$, $R_{t-1}$, $R_{t-2}$, ..., $R_{t-12}$  \n",
    "<font size=4 face=\"宋体\">默认参数：</font>\n",
    "    loss='ls',\n",
    "    learning_rate=0.1,\n",
    "    n_estimators=100,\n",
    "    subsample=1.0,\n",
    "    criterion='friedman_mse',\n",
    "    min_samples_split=2,\n",
    "    min_samples_leaf=1,\n",
    "    min_weight_fraction_leaf=0.0,\n",
    "    max_depth=3,\n",
    "    min_impurity_decrease=0.0,\n",
    "    min_impurity_split=None,\n",
    "    init=None,\n",
    "    random_state=None,\n",
    "    max_features=None,\n",
    "    alpha=0.9,\n",
    "    verbose=0,\n",
    "    max_leaf_nodes=None,\n",
    "    warm_start=False,\n",
    "    presort='auto',\n",
    "    validation_fraction=0.1,\n",
    "    n_iter_no_change=None,\n",
    "    tol=0.0001,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# micQ1 = pd.read_excel('result.xls',index_col=1,header=None)\n",
    "# reanalysis = micQ1.index\n",
    "reanalysis = ['stl3', 'd2m', 'tcwv', 'tcw', 'stl2', 'mn2t', 'tsn', 'stl4', 'stl1',\n",
    "       'ro', 'swvl1', 'swvl2', 'swvl3', 't2m', 'swvl4', 'mx2t', 'sf', 'cp',\n",
    "       'tp', 'rsn', 'lsp', 'sd', 'smlt', 'istl1', 'istl3', 'istl2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "预见期1:  ; 最优降雨滞后：6 ; 最优再分析滞后：8\n",
      "Wall time: 46.4 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAADKCAYAAAA2EILZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAASkUlEQVR4nO3de3Bc5XnH8e+DjS1CQLZ8QbZkiDFkjGUzk0RJMORSWuIy4VLSNA4t1AkQOyGXFiadTJtMp+mFMENIISHhYmMuAgqdScqQQKdjaCgTIE0qJ+AbbYItsKS1hGTZxgQ2WObpH2etrOR3V2elc/am32eGsVbn7HMeX/TjPXvOe15zd0REZLRjKt2AiEg1UjiKiAQoHEVEAhSOIiIBCkcRkQCFo4hIgMJRRCQg9XA0s+lmdpOZLcm9vs7Mzst9faqZPZh2DyIipZqe9gHcfdjMDgBmZn8AvET0ohV4G3A49D4zWwesAzj++OPfs3Tp0rRbFZEpZvPmzYPuPi+0LfVwzNMAvB9YDpwOtAEGLDWzNnffnr+zu68H1gO0t7d7Z2dnGVsVkanAzF4utC31cDSzacBJwGx3/4aZfRrocfcnctvbxwajiEilleO0+jBwdd7re8ZsvzztHkRESqWr1SIiAQpHEZGAcl6QScWmTZvo6+sbeb1z506WLFky8rq5uZlVq1ZVojURqWE1FY4PrN9Admh/0X0O73+VQ5n+kdfdmX42/uL5kdcNTbO4bN3a1HoUkfpQU+H4wvYdtPq0ovu0Mm1UOI7VtSeTdFsiUodqKhzPaFs27sixa0+GxQsWFq7RNCvptkSkDtVUOIZOh8d+5jjtjdc4duFJI6/1maOITERNhWOIgk9E0qBbeUREAhSOIiIBCkcRkQCFo4hIgMJRRCRA4SgiEqBwFBEJUDiKiAQoHEVEAhSOIiIBFVua1czeZWZ/Z2YPm1lD2n2IiJQi9XB092HgqKVZgW3u/vfA80SrEI5iZuvMrNPMOgcGBtJuU0RklHKeVh9ZmvVcYJW7HzKzxcAP3f2NsTu7+3p3b3f39nnzgsvKioikpmJLs5rZKuCPgT1mlnX3HWn3IiISV6WXZt2U9vFFRCZCV6tFRAIUjiIiAQpHEZEAhaOISIDCUUQkQOEoIhKgcBQRCVA4iogEKBxFRAIUjiIiAQpHEZEAhaOISIDCUUQkQOEoIhKgcBQRCVA4iogEKBxFRAIUjiIiAQpHEZGASq5b/Q4zu9TMvmBmMwLv09KsIlIxlVy3ei3RAlt7gFWB92lpVhGpmNRXH8xzZN3q5cDpRMF8CDgMHLVutYhIJVVs3WrgZeCjwInAf6Tdh4hIKSq9bvWv0z6+iMhE6Gq1iEiAwlFEJEDhKCISoHAUEQlQOIqIBCgcRUQCFI4iIgEKRxGRAIWjiEiAwlFEJEDhKCISoHAUEQlQOIqIBCgcRUQCFI4iIgGxnudoZucBlwDHA3uBR939v1LsS0SkosYNRzM7Gdji7k/kfe8UM2ty96FUuxMRqZBxw9HddwOY2UpgEDjD3X+YdmMiIpVUymeOjcBfU+JiWPlLs5rZKjPbaGbvNbOlZvZ5M7vMzOaU1LWISMpKCcfZwEagrZQD5C/NCjwJbAZ6gYPAQqAL2D/2fVq3WkQqadxwzI32PgI86+7PuvvNkzjeYaJwXAvMBP4VuABYNHZHrVstIpUU5zPHTQC50+LzgLfc/cdxD5C3NOsC4FzgJaKlWN8OvA/YiW4pEpEqUzQczWxN7ssTgd+6+4ZSDzBmadafjNm8pdR6IiLlMN6I7SRgqbt/l+jzQRGRKWG80+qn3P3nua9/kXYzIiLVoujIMS8YcfcfAZhZU9pNiYhUWqzpgwBm9iCQBeYCF6XWkYhIFYgdjkAHsANoTqkXEZGqUUo4Xgm8AMwHfpZOOyIi1aGUcHwa2EZ0Wi0TsGnTJvr6+kZe79y5kyVLloy8bm5uZtWqVZVoTUTGiPNUnmPc/a3cy1aikePY70sMY4Ovo6ODNWvWFNhbRCopzsixxcxOdfdvA5jZsWb2IaJbe15Ltbsadvf6mzm4r7foPt2ZIb6zZ2vB7SfMbuGKddck3ZqIxBBn+mC3mc00s28Bc4geGnGDuysYi9ixfSsnWPFwPMFg/57C+3Rn9LhMkUqJ9Zmju78IfDnlXurKsrYVHNxX/JbQ7swQixYW3mfR7Jak2xKRmEq5ICMliHM6rM8cRapX3DVkVgM9wIvAx939tlS7qlOhq9UdHR0jr3W1WqR6xB05ZoAT3f0VM9uWZkP1TMEnUjvihuMFwGEzWwgs4+hHj4mI1JW44bgD2J37+nBKvYiIVI24V6vvO/K1mS1Nrx0RkeoQe3kCM7vYzL4PnJJiP1Unm83y0EMPkc1mK92KiJRRrHA0s9uInsbzLXf/aikHKLQ0a27b+Wb2HjNrKL318ujs7OSuu+6is7Oz0q2ISBnFPa2+2sw+ALSZ2UDupvBY3H3YzPKXZj0N6DWzs4GVwD6ip/2MYmbrgHUAJ598ctzDJSL/lpvh4WEWLVrEr371K3bt2jWyj267Ealv5u6lvcHs6lLvczSzrwP3A7uA9wLnA68CjwPvAl539x8Uen97e7unOXK7/45beGNoT8HtL+0Z5B0Lij+M6LimBVz+2S8l3ZqIpMjMNrt7e2hb3JvA/4FoSVWIFt0q5eCFlmbtBy4EhnOvK+bFl7rZv7e/4Pb9B1/nN9nfFq0x69XhpNsSkQqKeytPlijcni71AOMszfrdUuul4evX3zDq9XjPXQSdVovUu7ifOX7DzJqBDwC/TLelylPoiUjsW3mI1q3+BPDBlHoREakacT9zXA+8AnyO6JmOIiJ1Le5njj8lutJ8JnAq0dN5RETqVtzPHO/Oe/lUSr2IiFSNUj5zFBGZMhSOIiIBCkcRkQCFo4hIgMJRRCRA4SgiEqBwFBEJUDiKiAQoHEVEAhSOIiIBCkcRkQCFo4hIgMJRRCQg7iPLZAoZu0wEHL1UhJaJkHqXejia2XTgm0TrxSwBPgnc7u7/Y2YLcl//UeB9FVuadaoLhV5HRwdr1qypQDcilZF6OBZZt3oOsAg4WOB964H1EC3NmnafU9mtd97EwIHeovv09w7RNbCl4PZ5jS18/jPXJt2axDDegnATGeXr7KH8p9WHgc3AWmA30AgsNbOV7v7TMvciOdt2bOX1GT3Fd5oBXQOF9+nvHUq4K4lrbEAlMcrX2UN5TquD61a7+89y29sVjJW1fNkKBg40Fd2nv3eIk1oK7zOvsWXc46QxwhFJSzlOq4utW427X552D/UsiTW245wOpzEamWojEaktulpd49IInLGBm81mefzxx8lmszQ0NADxRnk3b7yVngOvFNw+1PMKWwZ3Fa3R2jifa676fAnd17c7Nt7P4IFs0X36erroHryz6D5zGxv47FW/G5fcveEBDu4vXrc708V3+jcW3H7CrAauWHtZ0RppSGKAEKJwrDEbNtzEvn2FP/vLZPbR1/d80RqzZ7eydm3h0eLYf0QbNmzgwIEDmFlJwbtl+zZ6G35TeIcG6B3sKlpjqKdwuE5FgweyHLf4wqL7LF4co07Xo6NeH9yf5cwFFxd9z5kLitfcsueH4x84BWmdkSgca8z27Vs5ZlrhcDxmGvT1F7+4ksnsK+mYq1evZteuXaxevbqk9zW1zqfv1DdLes9RNXbNmNT7601fTxfDg/826TrTs6PvTujOdLG//weTqnnw8DgX9WqMwrHGtLWtYN++2QW3ZzL7WLiw8HaIRo6laGxs5Prrry/pPRCdErOr+Gl1U+v88WvIiOVtZ8Q6rW5uLT58nNt4xqjXy5adEeu0etHCwnUXzTqj4LYkPXDHPbyxN3gHIABde7q5s+eWojWOm3MCl33200X3Mffqv4Wwvb3dOzs7K91GTailixy11GstSevPtVr+vr72pa/Q+taJk6rRc8yrXHfLDZjZZndvD+2jkaOI1JRly5eNO3JcvGBR8Rpzxj97UjjWuNCVuo6OjlH7VMv9g+P1Wi19SnUb73Q4qRGuTqtFatxUnz6Y3+vw8DDPPPMM55xzDtOn/27sV6hXnVaL1LE0AqoaQi+u/F6ffPJJHnroIRYvXsy55547qbp6nqOI1I1Dhw6N+nUyFI4iUjeOPfbYUb9Ohk6rRaRs0n74yMqVK7nyyitZuXLlpPoEXZARkQqq9L2TuiAjImV3/+33kx0sPuumq6+LO3cXf0hGw9wGLv9c+R/epXAUKZNauj0mCS9se4GWQ8Wf89lCC292F59/39VX/OEkaVE4ipTJVHu69uLmxVww44JJ13nszccS6KZ0CkeRlHx7471k9hee5gawt3c32we+W3D7wlkn8JdXfSrp1sqiYW4Djw0WD7auvi4WNxd/SEbD3IYk24pN4SiSkq3bt9M3s/gTkpg5m77BwgG6t3d3wl2VT5zPCat55FyxpVmBfcAngJXAWnfvT7sXkXJa0dbGnBgjxzkthZceXjhrYdJtVVQtza+v2NKswIC7X29mXwSOeqKp1q2WWhfndLiaR05pqJbgi6PcM2RGlmZ190Nm1gz80t27x+7o7uvdvd3d2+fNm1fmNkVkqqvY0qxm9m7gC8BLZoa7P5N2LyJxlfNJN9V6WjnVaYaMSAxT7fR3qtAMGalbU+3GaikfhaPUtKRurL75zo30HjhQcPtQbw9bBwaL1mhpbOSaz1xV0nGleikcpabcvHE9vfuLLy071Jth68CegttbZs3mmqvWjfre1h3bycwocrPxjAYy44TjUG99LU061SkcpaZs3b6d3plWfKeZRu9g4XAc6s0c9b0Vy9poGmfk2NRSfFGmlsbG4n1JTVE4Sk1Z0dZGU4yRY1NL4ZunW2YdPWtlvNNhXZCZenS1WupOEkGW9kNZpTroarVMGdlslm3btpHNZmlomPgDCxR8onCUmjZ2hLd7926ee+45brzxxpFppxrlyUQoHKWmjQ29bDbLaaedxiWXXDKpkaOIwlHqSkNDA5deemml25A6oKVZRUQCFI4iIgEKRxGRAIWjiEiAwlFEJEDhKCISoHAUEQlQOIqIBKQejmY23cxuMrMlZrbKzDaa2XvNbLaZfdrM1ppZU9p9iIiUopJLs/4J0UqEAH9GtK71iPylWYHXzOz/Yh5yLlD8qaSlS6NmWnXVq3qd6r2WUvOUQhvKPX1wZGlWYC9wiGj0+sbYHd19PbC+1AOYWWehRxBNVBo106qrXtXrVO81qZoVW5oV2AlclNvt4bT7EBEpRTlOqw8DV+de/mTM5rvTPr6IyETU49Xqkk/FK1QzrbrqVb1O9V4TqVkTyySIiJRbPY4cRUQmTeEoIhKgcBQRCaircMyfjZNgzQ+b2XVmdl+CNU8zs2vM7J+SqplXe4WZbUyw3nG53//fmpklWHeama02s7YEa37MzK43s0fMrDXBun9qZh8xs+KLW5dW8xNm9odm9rWE6uXPRLvczM41swsTrjs3iZ+DMTW/aGZ3m9llCdcdmY030Xp1FY7uPgwcmY2TlGfd/WtAJqmC7v4i8GPgf5OqCWBm7wR+A8xMsOyZRPepvujJXr37K2AacHZSBd39YXf/G+ARd+9Jqi7RhIVTgOcTrPkx4D+BE82sebLF8v7tvw14t7s/CXw8wbrm7oMk8LOVV/MY4Hbgy0BLgnWPzMbbTDQbb0LqKhzT4O6HzGw5cFfCpfuApNcL/SRwMbDUzBYlUdDdf+bunyHBEMt5p7s/CJxlZnOSKpoLmsT+R5bnUeAKM0vqZ+Zmor+v04GBhGpC9DN9OPf1awnWTYMT9Xo2cFPCtfNn401IXa0+OGY2zosJ1bwceB+w38y+5+79CdQ8D5gBPDLZWvnc/R9z9dvdvTuJmrlTswaicEjSw2Z2EbDV3fcmWPcjwL8kWA9gPrCMaHZXUjLAQuCruYkSk5L3b/9EYIeZfQh4LMG6C8xsLzDfzN7u7hMO3jE/p9cC/UQj8+8l1SujZ+NNrJ7ucxQROZpOq0VEAhSOIiIBCkcRkQCFo4hIgMJRRCRA4SgVlZuBdI+Zxb6tzMzebWYfLrDt2hLqXJzkbCqpL3V1n6PUpJeBHmCemd0C/ILoKfHvB7YQzST5C+A64Fbgo8Druf+eAjCzK4A3geeA883sUeAGovv8XieaNTIf+DnQBLxFdIP0hcA0M/sK8M9ENyJ/CrgX+FHu+KcBw+7+QJp/CFJ9NHKUquDue4As0Q3cy4FXgaeBJ4B5RLMpjsx8yQDb8t6+F/i93K8D7v5romlk9/C76WNH3numu3+faHpZBtia+/Uw8AqwDxhy99uAPwfuA/490d+s1ASFo1TaqUBLbrpjE9BMNMsB4Cyi2TlbgOOBDwOziaZenpX3IIwVwOPAycDruYdZzAUWEc2S+H2iaXpzgGPN7OtAK1EongPsANbkjtECtJrZ8cBPgQeAD6bzW5dqphkyUnVyoXcrcK+7/3el+5GpSeEoIhKg02oRkQCFo4hIgMJRRCRA4SgiEqBwFBEJ+H/4OY++YfdzFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAADKCAYAAABjXaPCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdgklEQVR4nO3df3Rc5X3n8fcXnKxcY0saG9vSyD8EDocEzO4GJZt1nLJkE9JsCHWbU8ICJcG/gJ5uCyc9m2ybzaa7pOFsTotpCD+MbYiDC4cNweFHyabZssEc0rAW3djgpFnM2JZmJFvSSIptLDuWn/1jJHsk3TvS3LlXc+/o8zqHI808er7zlbie79znPs9zzTmHiIiIJMN51U5AREREpk6FW0REJEFUuEVERBJEhVtERCRBVLhFREQSRIVbREQkQVS4RUREEiTywm1ms8zsXjO72MyuMbOtZvaBkbY/MrN7zGxx1HmIiIjUgllRv4Bz7rSZDQIGvASsALJmlgIuAQaBw+P7mdlGYCPAnDlzrrz00kujTlVERCQ22tvbe51zF45/PvLCPc4w0A5scM79OfCHZrYeWAnsKf5B59xmYDNAW1ub27179zSnKiIiUj1mdtDr+cgLt5mdDywCmoCrgQPAD8xsJXAlcBR4M+o8REREasF0DJUPA3eMPNw1rnlv1K8vIiJSSzSrXEREJEFUuEVERBJEhVtERCRBVLhFREQSRIVbREQkQVS4RUREEkSFW0REJEFUuEVERBJEhVtERCRBVLhFREQSRIVbREQkQVS4RUREEkSFW0REJEFUuEVERBJEhVtERCRBVLhFREQSRIVbREQkQVS4RUREEkSFW0REJEFUuEVERBJEhVtERCRBVLhFREQSRIVbREQkQVS4RUREEkSFW0REJEFUuEVERBJkVrUTEJGZ6Yc//CHd3d0A7N+/n4svvvhs2+LFi7nmmmuqlZpIrKlwi0hVFBfm7du3c8stt1QxG5Hk0FC5iIhIgqhwi4iIJIgKt4iISIKocIuIiCSICreIiEiCqHCLiIgkiAq3iIhIgmgdt4hMi01btpIdHPRsy2c72dvT69mWrq/nzvXrokxNJFEiL9xmNgv4BnA/cDHwWeAh4JfAZ4BrnHM3ePTbCGwEWLp0adRpikjEsoOD/GJZq3fjslaO+HU8mIkqJZkGQ0ND7Ny5kzVr1lBXV1ftdGpC5IXbOXfazAYBA14CVgBZ59ygmT0PNPr02wxsBmhra3NR5ykiM5O2Xo3W7t272bZtGy0tLaxevbra6dSE6R4qHwbagQ3AnwPHgQVmdolz7pfTnIuITKN8tpPTPsPhJfudGoogm3O09Wq02traWLt2LW1tbdVOpWZMx1D5+cAioAm4GjgA/MDMLgH+BfAG0BF1HiJSXal0C0f8hspL9dNQeaLV1dVxww0TroZKBaZjqHwYuGPk4a5xzTrLFhERKYOWg4mIhGhoaIgnn3ySoaFoh/gluSo9RrQcTESmRbq+3neGeD7bSSrd4t/PRxwnlmkyVm0LY5Z8pceICreITItSa7GDTgqL48SylStXcuWVV7Jy5cpqpxILcVoOFoeiC5VP2FPhFpEZ576t3yY3cHTC833ZQ7zZc79nn+aGufzxus9NGru9vf3sfx/96EcrzrWaBgcHueeee/jSl75EfYmRj1LiNALx8ssvs23bNlKpVODRmDBmyVc6YU+FW0RmnNzAUd5a/sGJDcs/SL9fpwOv+cYrHrLPZAqXA3bt2kVnZ2ei14I/9dRTtLe389RTT7Fhw4ZAMeK0HGz//v1jvgYRh1nyKtwiIhUqLsxDQ0N85Stf4Ytf/GLVh4Yrdf311/P2229z/fXXB44Rh0I36sYbb+TQoUPceOONVXl9vzkZ5X64U+EWkRmnL3uIk70Th8pL9jnpey4+Rl1dHZdffnnVi3YY13Pr6+v5+te/HnJm1VPt3yesORkq3CJSFePPPrZv3362Lerh5fnppfR7DZWX6lNiqDxscZlEFZeJZWFca69E8bEKlZ0th0GFW0SqIsib3aYtj5Ib/JVnW1+2gzd67vNsa66fx53rby379aolDjOXAV599VW2bdvGwoULqzrR7oknnqC9vZ0nnniC22+/PXCcoB9Exh+r1V7BoMItIomRG/wVv1z2Pu/GZe+jz6/jwX1RpRSJMJaUhXFt2czGfC1HmGvsW1tbx3wNKk4z3Cuhwi0iM05zw1zPWeJ92UPMT3vfRri5YW7UaZ31+uuv097ezuuvv87VV189ba87arTonj59miVLlpDJZOjoKNxSYqpFN8w19ldddRX9/f1cddVVgWNAvGa4V0KFW0QSoy/bwameqU0SG9Pv1LExj/3WY5dTYO7b+h26Bk54v172APt6Nk94vqlhNn+87vc9+3gtKXv55Zfp6Ogo6ww1jJnLxT83a9asqg0Lj/9dTp06dbYtyLXlOM1wr4QKt4gkxvz0Evr8hspL9YtgqLxr4AQHWj/u3dgKnnPWM3/nGy+sJWVx3E0uqDB+l7hNLAuDCreISMzEZUlZLYjbxLIwqHCLSNnC2kiiXM3183wnmvVlO5ifXuLfL2R92QO80/tceX1Odk/+QxILOx5+jBN93mv9M10dbOn8pmfb7Plzuem2z0eYmQq3iARQreHYUku6pvtMan56OUf9hsr9+pQYKg9q69YdDAx43x4ym83Q07PVs62hoY51624KPZ9acaLvKJ++4F96N77H53ngub5/jCijc1S4RUQCaGqY7XvNui97gPnp5d59QjYwMMTy5b/t2bZ8YgpnHTjw/dBzkemhwi0iEoDf7HBI3nXURx/ZwVGfs/aOXIa/Pux91j63oY5bN8TvrP3xh7cz1Hvcsy3TfZAtHQ96ttUtmMPNt8X//5sKt4hIlTy89XF6B70LZndnho7eLZ5tC+rruG3dzaHlcXRgiCuarvNsu6LJv9+ermfPfr998w6O572Xxx3syvBgzvt3mZOazS0bwy3+Q73HuXb2R7wbW32eB57v3XX2+0xXB987M1j2a3ee572zX5hUuEVEqqR3cIjZrdd6tpXaJKw38/zZ77PZDL293y37tU+e7Cy7TynH8ydYnfIu/qtT/v1eyT875vHjDz3OUO/EDzOZ7gxbDnkX/7oFddx8e3gfZABam5b4X+Mu4bljusYtIhKpat7sJAyXXfbekpPT0mnvTwANDe+NMq3AhnqH+NS7PzWxwXtDOwBe6H1hzONM90G+d9p7qLyU7Kzes9/Pnj/Xd6JZpquD1ibvFQyz55/bYW/Hw9s4kfc+a890dbIle693jFQ9N9221jfPKRVuM/sYsAaYA/QBzzvn/vdU+oqIxFncC/NkSs0Mn+q19o5choHDT5f92keHz521H+zKcLyr/Bi9hHvmD9C6eJn/UHkJz584N1ReaknXVP+uJ/KDXDf3Eu9Gv+eBZ/O/LBl30sJtZkuBPc65HxU9t8zMUs65/GT9RSTcGy5Uw6YtW8gOem81ms/m2Ntz2LMtXd/InevXR5bX6N/12LFjvPjiixw7dowLLrgAmN6/q9fuXKNn7qXy6O7McLr3e2W/3qyhbLBEfSxpbvW9xl1K8TXuZU2tvkPlpYwfKg9D3YI5Y65XF8t0H6R18TLffkkwaeF2zh0CMLN/DfQC73XOhf+XFqlhSd+GMjvYz8+XNXs3LmvGu2wDB3NRpQSc+7t++ctf5uTJk+RyOe6+++5AsSq593TQDwiXX/bekpPTFrd4D3MvqA93mHtuQ92YIlysI5dhSbN3HnMbwt/ZLdOd4ZlfP1NWn+y7xn6QKTUzPIn//sYr5xp3PbAeeDKiXEQkpvLZHKd7jpTf79TpCLKZqK2tjddee62iuz5V45aPpWaGT2eBKbWka6p5zEnN9j17PtiVYVmTd/Gfkxq7tr11cav3Ne4SXjj1wuQ/VEPKKdyNwFbgg0D42/+ISGyl0s0c9jvjLtUv4jPuUfPmzRvzNYhaueWj32S7qC8dlFrSVQtnuUFkujp5Olv+FeXO894p2T6Va9zXAA541Tl3EHi17CxERCK0atUq1q5dy6pVqwLHqJVbPsZ9vsRk6hbUTZglDoUh9NbF3mftdQvieTOW1qYW/8lpJTx7tMLJac65HwKY2cUjs8vPOOf+vuxMRCSx0vWNvter89kcqbT32Xi6vjHKtM4Ko+hWco07DEEnuEWdy3QvkfNbj13Ns/Y4/b+BSQq3mY3+leYBJ51zj0SfkojETamZ4bUyDFqNa9zF4nSmHKdcKhVG0Y3b32OyM+5FQKNz7k/N7L9MR0IiItVQK9e44yIuG9vEreiGYbLC/WPn3Gsj378edTIiItUSp2vc1R62D0MtFsxyzU7V+26mkunqpLWpxbdfKSULd1HRxjn3HIA2XhERiVa1h+0lHKW2La3kEtOUl4OZ2RPAELAA+HSgVxORmlCtJUczhYbtpZRy1nFvB/YBiyPKRaSmbNr6AJ2DEzctyXceYU/v2559WuoXcue6P4g6tYqpOEcrTsP2Ej/lFO61wM+BhcBPo0lHpHZ0Dh5h30WnJjZc1EA3Hs8DvF3+7mQitSouE9zippzC/QrwBoWhcpEZI+k3CBFJKv3b8jaVndPOc86dGXnYQuGMe/zzpfrPAr4B3A9cDHwWeMg593/M7LPAbOfcYx79NgIbAZYuLXEjVpGIJf0GISJSW6Zyxp02s4ucc/cBmNm7zOw3KSwPOzZZZ+fcaTMbBAx4CVgBjN7K5fDIY69+m4HNAG1tbW4KeYqIiNS88yb7AedcB5A1s780s8eArwI/c85NWrQ9DAPtwIYAfUVERGa8KV3jds69BXwhyAuY2fkUdmBrAq4GDgA/GGluHnleRESkpoW1jLKcyWmBOOeGgTtGHu4a1/Y3Ub++SLnu3baJ7ECnZ1u+M8+e3j2ebemGFu5ae2eUqYlIgoU12W5KhdvMrgc6gbeAzzjnHgzl1UViKDvQyRsr9nk3roAc3d5tb419mO88wnDv8bJeOz80p6yfF5GZZ6pn3DlgnnPuiJm9EWVCIrUi1bKQbq913KX6vP3uiLIRkVox1cL9KWDYzJqB9zFuyFtERESmx1QL9z7g0Mj3wxHlIiIR87o38eiGMtpMRiQZpjqr/Duj35vZpdGlIyJRGl+YtaGMSPJMuo57lJldZ2bfBZZFmI+IiIiUMNVZ5Q8C/wj8pXPuJ9GmJFIbWuoXet40JN95hFTLQv8+IiIlTHWo/A4zWw1cZmY9IxuyiEgJfrfn1PC0iFRiykPlzrlXnHNbgI9HmI+IiIiUMNWh8v9KYatSKGxfKiIJsWnrZrID/Z5t+WyOvT1dnm3phkbuXLcxytREJICpLgcborCn+CsR5iISC/nOPGf6yl/1mD+RjyCbymUH+tm3POXduDzltw8cHIjn7yMy0031GvdfmNliYDWFSWoiNSvVkiK3wrec+fd7y6c4ioiEaMrXuCncyev3gI9ElIuIiIhMYqrXuDcDR4DbgfmRZiSJp925RESiM9Vr3D8B3gauAC5iwn2QRM7R7lzxks/mGO71noBWst9JF0E2IlKpqV7jfrTo4Y8jykUkFtINLb4fTfOdeVIt3tey0w0tEWYVXCrdTLff5LRS/TQ5TSSWpnrGLTJj3LX2Tt82jR6ISLWVMzlNREREqkyFW0REJEFUuEVERBJE17jF19DQEDt37mTNmjXU1dVVO52qKV7etn//frZv3362TcvbRGS6qXDHxPjiMLruGapXHHbv3s22bdtoaWlh9erV0/76cRHG376axT/d0Oi7fWk+myOVbvbvJyKxk9jCXWubfBTnG5eZy21tbaxdu5a2trZqp5J41TweS90oJC7HmohMXWILtzb5iF5dXR033HBDtdMQEZEiiS3cEo1aG8kQEak1M7pwq0hNFHQk45FH7qW/v9OzLZfrp7v7Z55tjY0tbNhwV/mJiojMUDO6cGu4PTz9/Z0sXLTXs23hIgDvon7kcHQ5iYjUIq3jDtHQ0BBPPvkkQ0ND1U5FRERq1Iw+4w6blk9JlHRpR0QgYYV7x+ZHGMoPeLZlunJszXmPu9alGrhp44aRGA9zIt/rE6ObLbkOz7bZqQXctPG2kvkldfnU1kc2MZj3HsruzOXp7doz4fn6VAvrNvjfjEPCF8alHa/iP7qmXMVfJBkSVbiH8gNcV9/k3ej3PPBs/ty9iE/ke7muYbb3Dza0lojhXeyLJXX51GC+k9YL3/Rsa70QYOK9nDM90eYk0VBhFkm+RBXuWvLtzZs4ns96th3syvNAznui15xUms9tDPdMtzOXp7f7TFl9hs7oXs0iItUw4wp3pqubp3Ony+7XaeH+qY7ns6yq/4Vn26p6gCOeba9GUC9bmlO0XjjxrLqUTE9qzONcrp/uADPEzwz3l99JRGQGS1ThznTleNrnOnYpnTZ89vvWpsX+Q+UlPDtwwrctjvuMT7fm5kYWLvK+Tl7KkcPaD7uUTVsfJjvgfZkmn+1mb8/EORnphgXcua70fAwRSa5EFe7Wpmb/a9wlPDtY3tlkueK4z3g56lMtvtesO3N5WppTE56vT7WMedzY2OK7JjuX66e52btANza2eD4vBdmBXva1+nzQbG2l2+v5zOTzMUQkuRJVuOtSDWMmmhXLdOVobfK+y1FdquHs97NTC3wnmmW6umltWuzZNju1oMxsSzvYledYgCH7Pgt/rLzU7PCpfhAptftZEj/MiIjEVeSF28xmAd8A7gc+CVwJ/AhoB1YBTcADzrn+cf02AhsBli5dCnB2SZeXqRaHUku6prPALGtKsare+zp2Ka8OTjz7ldqVz3Yz3FveB7z8yUR9HheRMkX+L9w5d9rMBins0vYQMA9YT6Fgv0xhrdEngCfH9dsMbAZoa2tzUecpBVrnGy+p9GK6/YbK/fpk/OdjiEjyTedHcwcMUzjLvhdoBK4HVgJ/NY151KwwJsmpMIuIxNt0DJWfDyyicIZ9F3AYWAbsAN4Efu6c+6dy4+rMcKKkT5KTidINC3wnm+Wz3aTSE+dkpBvCnY8hIvEyHUPlw8AdIw93jWt+KWjcahbmxx/+Jid8Jskd6Orlkaz3+uzZqSZuvu0/AIWNVPzWZB/syrOsyfta9pxUuvyEJbFKLevShzORmWlGz2IJetZ+It/FJ+d6LsSBuYD3Ih1eLCrUpXY/0xtyeLTGXkRqzYwu3HrTrn26fCAitWZGF+6gDnT1sjM7WHa/rvPC/3M/unkTR/sn7nnekcvz113e+53PbUxza8j7nUv0NK9DRECFO5DlTQv45NzyN0958Wj4k4aO9mdZOX/fhOdXzge/Ifu9faGnIdNAhVlEQIVbIqAzQxGR6KhwS+hUmEVEoqPCHcDsVNOYGeLFDnT1srzJe0h8dqr8G6RMpiOXZ6BrePIfLHLU6V7aIiJJpcIdwOhabC/TPXN5SXOKlfN9lqb52Nun/c5FRJJKhTvh5jamPSebdeTyLPG4Hedon7jzuk4+uga71HXyB7bcS8/gxFn2AIezeTI9ezzbLqxP8wfr/e9wJiISFyrcCee3rCvpa5bHF+ap/j49g1lOLX3Ds61xKZwi593vUPk5iohUgwp3TIzf4Wt0FjZoJraIiJyjwh2CMIquCnM4DmfzvNNzpux+v3FKE/ZEJBlUuEOgohsfi9IpTi31Hg4v5d2HNGFPRJJBhbuGaLhdRKT2qXDXEBXmwuxwv4lmh7N5FqW9z6wvrI//THsREVDhDtXQ0BA7d+5kzZo11NXVVTudGanUkq6kz7QXEQE4r9oJ1JLdu3ezbds2du/eXe1URESkRqlwh6itrY21a9fS1tZW7VRERKRGaag8RHV1ddxwww3VTiOxHt5yL30e9xYH6Mrl6TzsvevZ/MY0t2nXMxGZIVS4JTb6+rPMSe/1bFuRBvAu6n3eTwOaaS8itUeFW2qaCrOI1Bpd4xYREUkQnXFLbHTl8pw54srud95pbVcqIjOHCrfERlNzijnpEhesfRzPartSEZk5NFQuIiKSIDrjltiY35j2nSHelcvT1Ox9Zj2/UduVisjMocItsVFqLba2KxURKdBQuYiISIKocIuIiCSICreIiEiCqHCLiIgkiCanSSwV7zEOY/cZ1x7jIjKTqXBLLKkwi4h401C5iIhIgqhwi4iIJEjkQ+VmNgv4BnA/8EngSuBHwOvAeqAO2OWcezLqXERERJIu8sLtnDttZoMUzu4fAuYB651zO4AvmNlNwDPj+5nZRmDjyMNjZvZPk7zUAqC3wnQVI765KEZ8c1GM+OaiGPHNZSoxlnk+65yL/D/gq8AKwIBrgXcVtd0a0mvsVoxwY8QpF8WIby6KEd9cFCO+uVQSI/Jr3GZ2PrAIaKIwXP5+Rs6kzexDwKtR5yAiIlIrpmOofBi4Y+ThrnFt/xD164uIiNSSWppVvlkxQo8RVhzFCD9GWHEUI/wYYcVRjPBjhBWnqjFsZKxdREREEqCWzrhFRERqngq3iIhIgqhwi4iIJEhNFG4zm2Vm95rZxQH7X2VmXzOz71SYxwozu9PM7q4wzkoz21pB/9kjv89/NjMLGON8M7vezC6rII/fMbOvm9n3zaylgjj/3sw+bmbrKojxe2b2CTP7s4D9zx5jZnazmV1tZtdWEGNBkONtXIw/NLNHRzYxChrjGjPbamYfCBpj5PHXzOxjFeTxYTO7x8w+VU4MjzjvMbPfMrNFFcT4mpndbWaPVBDjy2b2ATP77Qpi/ImZfSZAjLPvZxUcq8Uxgh6rxTGCHqvFMYIeq2Pe3wMeq8V5VHKsFscJdKxCjRRu59xpYJDCBi9BvOqc+zMgV2EebwF/D/wiaAwzuwQ4DvyzClK5gsLa+bdc8NmHfwKcD6wKmoRz7hnn3H8Cvu+c6wwaB+ijsIPQzyqI8TvA/wLmmdnicjsXHWO/AbzfOfcS8JmAMcw510uA47UoxuhOhF8A0kHzAF4C2oFs0Bhm9m+BA+X098jj/cBCYF+Fcb4KzAeWBI0x8l5wH/D9CvJ498jXsn6fohiNwEXOuaeB3y0nBmPfzwIdq8Uxgh6r4/IIdKyOixHoWC2OEfRYHZdH4GN1XJyvEuBYhRop3JVyzv3azC4HtoUQrhuo5J6UnwWuAy41s7L/hwI4537qnFtPBUUXuMQ59wTwITObHzTISJGs6APRiOeBW80s6DG7icLf9j1ATwV5nAcMj3x/rII4lXIjeawC7q0gzjCFN8MNAfvXAf8KuJoKjnvn3DcpvLF/PmgMCsWuC/gfwBcriAPwCeB/VtB/N4UPeR8N2H8A+JmZ3QxkyulY9H72NwQ8VsN4TxwXI9Cx6hGj7GN13N8j0LFanEclx2pRnO9RwbFaE/fjtrG7s70VoP/NwAeBATP7lnPucMA8Pkbhk3ZZn9SLOef+20isNudcR8A8rqXwhvp80DyAZ8zs08Be51xfBXE+TuEfTCUWAu8j2CflUTmgGfjTkU2BylJ0jM0D9pnZbwIvBIzRZGZ9wEIzu8A5N+U31XHH+l3AYQqjEd8KGONqCn/XH0y1/7gYjc65vzCzzwNljaqM+3tcS+Hf7hPlxBgXpw44Any43Dge7yHnO+d+XUGMf05hhOdEBTF+SmH07Gtlxjj7fgb0BjxWx7wnEuxYHY1xEriAQqEq91gt/l2OU7g5VbnHanGMb1G42VW5x+pojEEzOwa8SbBjtTiXHgIcq6B13CIiIomioXIREZEEUeEWERFJEBVuERGRBFHhFhERSRAVbhERkQRR4RapESO7Mj1mZlNe5mlm7zezq3za7iojznUWcOdCESlPTazjFhEADlJYn3qhmX2TwprX/RQ2ndhDYfe4P6KwLvgB4N8B74z892MAM7sVOAX8X+C3zOx54L9TWAf8DoVNRRYCrwEp4AyFzT2uBc43s/8I/BWFjTY+B3wbeG7k9VcAp51zO6L8I4jUOp1xi9QY51wXMERh45vLgV8BrwA/Ai6ksOva6G52OeCNou59wL8Z+drjnPt/FLbgfIxzW02O9r3COfddCltR5oC9I1+HKWyE0g/knXMPAr8PfAf421B/WZEZSIVbpHZcBKRHtspNAYsp7MIF8CEKO4vtAeYAV1HYIrSbwra2o3tRrwT+DlgKvGOFm8wsoLCf8gEKW3i+h8Iey+8ys68CLRQK9ocp7N98y8hrpIEWM5sD/ATYAXwkml9dZObQzmkiNWykID8AfNs59w/VzkdEKqfCLSIikiAaKhcREUkQFW4REZEEUeEWERFJEBVuERGRBFHhFhERSZD/DwwNnYBMn5e4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%time\n",
    "'''\n",
    "针对每个预见期寻找最优的输入。\n",
    "'''\n",
    "rcParams.update({'font.size': 7})\n",
    "Rlag = 13\n",
    "itera = 50\n",
    "model = GradientBoostingRegressor()\n",
    "for leadtime in [1]:\n",
    "    df = pd.DataFrame(columns=['input structure','mae'])\n",
    "    for i in range(1,Rlag+1):\n",
    "        X, y =  MyPrepaData.Preparation(data=data_train,Qlag=[1,4],Rlag=range(1,i),isreanalysis=False,leadtime=leadtime)\n",
    "        length = X.shape[0]\n",
    "        tscv =  iter([(np.arange(length-731, dtype=int),np.arange(length-731,length, dtype=int)) for i in range(itera)])\n",
    "        score = cross_validate(model,X,y,cv=tscv,scoring=make_scorer(mae,greater_is_better=False),n_jobs=-1, verbose=0)\n",
    "        df = df.append([{'input structure':i,'mae':-_} for _ in score['test_score']],ignore_index=True)\n",
    "\n",
    "    df2 = pd.DataFrame(df.set_index(['input structure']).values.reshape(Rlag,itera)).median(axis=1)\n",
    "    opt_Rlag = df2.idxmin()\n",
    "\n",
    "    plt.figure(figsize=(5,3))\n",
    "    plt.subplot(111)\n",
    "    ax = sns.boxplot(x=\"input structure\", y=\"mae\", data=df,linewidth=0.5,fliersize=1)\n",
    "    plt.xticks(range(0,13),range(1,14))\n",
    "    plt.ylim(130,144)\n",
    "    plt.yticks(range(130,145,2))\n",
    "    plt.ylabel('MAE (m$^3$/s)')\n",
    "    plt.xlabel('Input structure')\n",
    "    plt.savefig('inputrainfall%d.tif'%leadtime,dpi=500,bbox_inches='tight')\n",
    "\n",
    "    Relag = len(reanalysis)\n",
    "    dfmic = pd.DataFrame(columns=['input structure','mae'])\n",
    "    for i in range(1,Relag+1):\n",
    "        X,y = MyPrepaData.Preparation(data=data_train,Qlag=[1,4],Rlag=range(1,opt_Rlag+1),\n",
    "                          reanalysis=reanalysis[:i],isreanalysis=True,leadtime=leadtime)\n",
    "        length = X.shape[0]\n",
    "        tscv =  iter([(np.arange(length-731, dtype=int),np.arange(length-731,length, dtype=int)) for i in range(itera)])\n",
    "        score = cross_validate(model,X,y,cv=tscv,scoring=make_scorer(mae,greater_is_better=False),n_jobs=-1, verbose=0)\n",
    "        dfmic = dfmic.append([{'input structure':i,'mae':-_} for _ in score['test_score']],ignore_index=True)\n",
    "    dfmic2 = pd.DataFrame(dfmic.set_index(['input structure']).values.reshape(Relag,itera)).median(axis=1)\n",
    "    opt_Relag = dfmic2.idxmin()+1\n",
    "\n",
    "    plt.figure(figsize=(8,3))\n",
    "    ax = sns.boxplot(x=\"input structure\", y=\"mae\", data=dfmic,linewidth=0.5,fliersize=1)\n",
    "    plt.xticks(range(0,26,1),range(1,27,1))\n",
    "    plt.ylim(127,135)\n",
    "    plt.yticks(range(127,136,2))\n",
    "    plt.ylabel('MAE (m$^3$/s)')\n",
    "    plt.xlabel('Input structure')\n",
    "    plt.savefig('inputreanalysis%d.tif'%leadtime,dpi=500,bbox_inches='tight')\n",
    "    \n",
    "    print('预见期%d: '%leadtime,'; 最优降雨滞后：%d'%opt_Rlag,'; 最优再分析滞后：%d'%opt_Relag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}