CEEMDAN_VMD_GRU_exmaple_of_sh000001.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "77c61528",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4/4 [==============================] - 1s 20ms/step\n",
      "======Co-IMF0 Predicting Finished======\n",
      "          r2      rmse       mae      mape\n",
      "0  0.939303  0.017272  0.013591  0.026389\n",
      "Running time: 107.557s\n",
      "4/4 [==============================] - 1s 17ms/step\n",
      "======Co-IMF1 Predicting Finished======\n",
      "          r2      rmse       mae      mape\n",
      "0  0.948441  0.010116  0.007522  0.018516\n",
      "Running time: 117.594s\n",
      "4/4 [==============================] - 1s 17ms/step\n",
      "======Co-IMF2 Predicting Finished======\n",
      "          r2      rmse       mae      mape\n",
      "0  0.977563  0.006268  0.005997  0.012192\n",
      "Running time: 118.332s\n",
      "======sh.000001 Predicting Finished======\n",
      "          r2       rmse        mae      mape\n",
      "0  0.975644  30.793454  26.876594  0.008086\n",
      "Total Running time: 463.494s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABZzklEQVR4nO2dd3iUVdbAf2cmvVcgjSSEQOgBIkVAEARBUewV+66rq+tW27or31p2dXWLrm2xomvviKAoghTpLUBoAdJ7b6TO/f5434QEEjKBFMjc3/PMk3nPLe+9M3DeO+eee44opdBoNBqNY2Dp6QFoNBqNpvvQSl+j0WgcCK30NRqNxoHQSl+j0WgcCK30NRqNxoHQSl+j0WgcCK30NQCIyP+JyP96ehy9BRFJEZELzPd/FJHXTrGfPSIyrTPH1lOIyDQRyejpcTg6WulrTgkRiReRrSJSZf6Nb1YmIvK0iBSar6dFRLqh7fkislJESkUkpZ3xR4mIEpEK85UiIg91zqfTEqXUX5VSP2uvnoi8JSJPHNd2mFJqVWePSURuFZEGc+5lIrJTROZ29n3aGUPTg1HTfWilr+kwIuICfAn8D/AHFgFfmnKAO4HLgFHASOAS4Bfd0LYSeAO4vwPT8VNKeQHXA4+KyOxW5uvUgf7OJtabc/cDXgI+EBG/Hh2RputRSumXA72AB4FMoBzYD8ww5f8HfAS8bZbtARLa6GOW2Yc0k6UBs833PwF3Niu7A9jQ1W2byS4AUtr5HKIABTg1k20G/gBMAzLMzyoHeAdjgfQQcAgoND+rgGZtbwJSzbJHgBTggmaf7f+a1Z1szrMESAduxXjY1QG1QAXwlVn3+H7a/I6AMcB2s+xj4EPgiTbmfyuwttm1h/l5nGNeuwLPmp9vLvAK4G6WBQFLzPEXAWsAi1mmgIHN+n2rcQyNn6v5/h3ABhw15/tAT//fcJSXXuk7ECIyGLgX4z+2N3AhhlJp5FLgA4yV32LghTa6GgYkKvN/r0miKW8s39msbOdxZV3V9pQwTUqTzH62m+J+QAAQiaGQf4XxC2QqEAoUAy+a7YcCL2Mo/lAgEAhv416RwDLgP0AwEA/sUEotBN4F/q6U8lJKXdLGcFv9jsxfO59jKNkA4H3gcjvnbwVuw3jopJrip4BB5vgGAmHAo2bZ7zEeisFAX+CPGMrebpRSN2E8UC4x5/v3jrTXnDpa6TsWDRgruKEi4qyUSlFKHWpWvlYptVQp1YCxEhvVRj9eQOlxslLAu43yUsDLtM13ZdtToQBjtfoa8JBSaoUptwELlFI1SqmjwF3AI0qpDKVUDcaq+yrT9HMVsEQptdos+7PZvjVuAL5XSr2vlKpTShUqpXZ0YLxtfUcTACfgebPfz4BN7fQ1QURKgGqMVf18pVSe+VnfCfxWKVWklCoH/gpcZ7arA0KASPNea457EGvOYLTSdyCUUsnAbzAUVp6IfCAioc2q5DR7XwW4tWHPrgB8jpP5YJgVWiv3ASpMxdCVbU+FIKWUv1JqiFLq+WbyfKVUdbPrSOBzESkxFeVejIdoX4zVfXpjRaVUJYaZpzUiMExEp0pb31EokHmc8k3n5GxQSvlh7I8sBqaY8mAMc8/WZvP9xpQDPAMkA8tF5HBXbYBrugat9B0MpdR7SqnJGEpMAU+fQjd7gJHNvWowNl33NCtv/ith1HFlXdW2Mzl+5ZoOzFFK+TV7uSmlMoFsDGUOgIh4YJh4WiMdiLHznh0hGwg77rOJaKtyi5sqVQHcDdwkIqMxfv0cBYY1m6uvMjZ9UUqVK6V+r5QagGFu+p2IzDC7q8J4YDTS72S3tmtmmk5FK30HQkQGi8h0EXHF+El/lLbNECdjFcYq9z4RcRWRe035D+bftzEUQZj5S+L3GLbmLm0rIhYRcQOcjUtxa+bZc7q8Ajxp2uQRkWARmWeWfQLMFZHJ5v0eo+3/W+8CF4jINSLiJCKBzdxOc4EBpzi+9Rifzb1mv/OAcfY2Vko1mrgeVUrZgFeBf4lIHwDz+7jQfD9XRAaaD5hS876N/452ADeIiNX0hJp6ktueznw1p4hW+o6FK8YGXQGGmaAP8LA9DUVkmYj8EUApVYuxqXkzhgfH7cBlphzgv8BXwC5gN/C1KevqtudhPMiWAv3N98vtmZ8dPIdhAlkuIuXABmC8Oa49wD3Aexgr7mKMjc4TUEqlARdhPMyKMJRk4y+b1zH2W0pE5IuODM78DK7A8HYqAeZjeNjUdKCbfwMXichIDM+lZGCDiJQB3wODzXqx5nUFxsPmJaXUSrPs1xhutiXAjcDJ5vE34E/mfP/QgXFqTgPR+y8aTe9ERDYCryil3uzpsWjOHPRKX6PpJYjIVBHpZ5p3bsHY7/imp8elObPorScNNRpHZDDG4S1P4DBwlVIqu2eHpDnT0OYdjUajcSC0eUej0WgciDPavBMUFKSioqJ6ehgajUZzVrF169YCpVRwa2XtKn3T73k1hrufE/CJUmqBiKzh2PH3PsAmpdRlpu/ucxhuaVXArUqpbWZftwB/Mts8oZRadLJ7R0VFsWXLlnYnqNFoNJpjiEhqW2X2rPRrgOlKqQoRcQbWisgypVTjkW1E5FOMkLcAczD8eGMx/JhfBsaLSACwAEjAOIm3VUQWK6WKT2VSGo1Go+k47dr0lUGFeelsvpp2f0XEB5jOsUMY84C3zXYbAD8RCcGI6PidGcCpGPgOOCF2uUaj0Wi6Drs2cs0j1TuAPAzFvbFZ8WXACqVUmXkdRstATxmmrC358fe6U0S2iMiW/Px8e+eh0Wg0GjuwayPXDOMab2bV+VxEhiuldpvF12PE7OgUzLjiCwESEhK0P6lGo6Guro6MjAyqq6vbr+xAuLm5ER4ejrOzs91tOuS9o5QqEZGVGGaZ3SIShBHUqXmyhkxaRvcLN2WZGJlzmstXdeT+Go3GMcnIyMDb25uoqChaBhJ1XJRSFBYWkpGRQXR0tN3t2jXvmNEE/cz37sBMYJ9Z3Jg8ovnjdzFws5mNaAJQap4K/BaYJSL+IuKPkfruW7tHqtFoHJbq6moCAwO1wm+GiBAYGNjhXz/2rPRDgEVmSjUL8JFSaolZdh1G1MbmLMVw10zGcNm8DYzQrSLyOEYeUoDHzHCuGo1G0y5a4Z/IqXwm7Sp9pVQiMLqNsmmtyBRGmNnW6r8BvNGxIXYdK/bmkl9eQ0JUADHBnvoflUaj6fU4bBiGIwWV3PW/rTz02S4u+OePjHn8O95cd6Snh6XRaByAVatWMXfuXAAWL17MU08dbzA5RklJCS+99FKn3dthlf5Ty/bibLXwyV0TefrKEQwI9uKpZfsorOhIzgmNRqM5RkNDQ4fbXHrppTz0UNtphrXStxObTdFga93jc8PhQr7dk8svp8WQEBXAtef05+krR1JTb2PR+jZPL2s0GgcmJSWFuLg4brzxRoYMGcJVV11FVVUVUVFRPPjgg4wZM4aPP/6Y5cuXM3HiRMaMGcPVV19NRYVxtvWbb74hLi6OMWPG8NlnnzX1+9Zbb3HvvUbmz9zcXC6//HJGjRrFqFGj+Omnn3jooYc4dOgQ8fHx3H///ac9jzM64Nqpkl5UxQ2vbeDhOUO4aERIizKbTfHE10mE+rrxsynH0nMO7OPFzKF9eXt9CndNHYCHS6/8aDSas56/fLWHpKyy9it2gKGhPiy4ZFi79fbv38/rr7/OpEmTuP3225tW4IGBgWzbto2CggKuuOIKvv/+ezw9PXn66af55z//yQMPPMDPf/5zfvjhBwYOHMi1117bav/33XcfU6dO5fPPP6ehoYGKigqeeuopdu/ezY4dOzplrr1ypR/q547NBu9uPHHV/tn2THZnlvHgnDjcnK0tyu6aOoCSqjo+2px+QjuNRqOJiIhg0qRJAMyfP5+1a9cCNCnxDRs2kJSUxKRJk4iPj2fRokWkpqayb98+oqOjiY2NRUSYP39+q/3/8MMP3H333QBYrVZ8fX07fQ69cjlrtQjXj4vg2eUHOFJQSXSQJwBHaxt45tt9xEf4cemo0BPajY0MICHSn9fWHmH+hEicrL3ymajRnNXYsyLvKo738Gu89vQ0dIxSipkzZ/L++++3qNdZq/TOoNdqtWsSInCyCO9vSmuSvfVTCrllNTxy8ZA23TN/MTWGjOKjLN2d011D1Wg0ZwlpaWmsX78egPfee4/Jkye3KJ8wYQLr1q0jOTkZgMrKSg4cOEBcXBwpKSkcOnQI4ISHQiMzZszg5ZdfBoxN4dLSUry9vSkvL++0OfRapd/Hx42ZQ/vy8ZZ0qusaKK2q4+VVycyI68M5UQFttpsR14eYYE9eW3O4G0er0WjOBgYPHsyLL77IkCFDKC4ubjLFNBIcHMxbb73F9ddfz8iRI5k4cSL79u3Dzc2NhQsXcvHFFzNmzBj69OnTav/PPfccK1euZMSIEYwdO5akpCQCAwOZNGkSw4cP75SN3DM6R25CQoI6nSQqaw8WMP/1jfz72nj25pSxcPVhlv16CnH9fIwKq58FixNMuBucXJvaLfophQWL97DkV5MZHtb5NjWNRtMx9u7dy5AhQ3p0DCkpKcydO5fdu3e3X7kbae2zEZGtSqmE1ur32pU+wLkxgUQFevDKj4d4a10Kl8eHHVP4h1bCD4/D9wvg5XONa5PLRofh6mRpYRo6HptNUVZd19VT0Gg0mk6ldyr9hno4sBxLRTY3jO/PvpxylILfzhxklNfXwrIHwD8arv8AlA3euQxWGafifN2duXhkCF/uyKKqtv6E7pOyypj34jom/e0H8spbBjvKK6/mue8PUtdg6+pZajSabiQqKuqMW+WfCr1T6ZdlwntXw/Z3uWpsBB4uVm6eGElEgIdRvvEVKDgAc56GwXPg7vUwaDZseMl4IAA3jOtPRU09S3ZmN3VbXdfAs9/u59IX1pJdepSqugZeWdXS9v/k13v51/cHWJdc0G3T1Wg0GnvpnUrfPxIiJ8PO9wnwcObH+8/n4YtMm1dZNvz4NAyaA4MuNGTObjD2NqguhSM/AjA20p+Bfbx4f7Nh4impquWGVzfwwspk5sWH8f3vpnL56DDe3ZhKXpmx2k/MKOHLHVkAWulrNJozkt6p9AHir4eiQ5C+iWBvV6wW00Xzuz9DQx3M/mvL+jHng6sPJH0BGP6314/rz/a0Elbuy+OqV9azO7OMF28Ywz+uGYWfhwu/mj6QepvipVWHUErx5Nd7CfR0YVSEH2uTC7t3vhqNRmMHvVfpD50Hzh6w871jsowtsOtjmHQfBAxoWd/J1TD17PvaeCgAV4wOw8XJwm1vbSa3tJpFt4/j4pHHwjpEBnpy5Zgw3tuUxrsb09h4pIjfXBDLrKF92ZtdRn65Dt6m0WjOLHqv0nf1hiGXwO7Poa4alILlfwLPPjDpN623GToPjhbDkdUA+Hu6cOWYcPp4u/LBLyYwMSbwhCb3nh+Lzab40xe7GRDsyXXj+jMlNgiAnw5pE49GozmzsCddopuIbBKRnSKyR0T+YspFRJ4UkQMisldE7msmf15EkkUkUUTGNOvrFhE5aL5u6bppmYy6HmpKYf9S2LcE0tbD+X8EV6/W68fMABevJhMPwOPzhvHTQ9MZFtq6v37/QA+uHBMOwMNzhuBstTAs1Bdfd2e77fptRQPVaDRnJkopbLaz00PPnpV+DTBdKTUKiAdmm7lvb8VIgB6nlBoCfGDWnwPEmq87gZcBRCQAWACMx0imvsDMldt1RJ8H3qGw7W34bgEEx8Hom9qu7+xmePHsXWK4fQJOVku7MXj+eNEQXrhhNBcMMU7ZWS3CuTGBrD1YQHuH34oraxn/1xXc/b+tlB7Vfv8azZlKSkoKgwcP5uabb2b48OHccccdJCQkMGzYMBYsWADA5s2bueKKKwD48ssvcXd3p7a2lurqagYMGHCy7rsNe9IlKqDCvHQ2Xwq4G7hBKWUz6+WZdeYBb5vtNoiIn4iEANOA7xrz4orId8BsoPUgFJ2BxQqjroW1/zKub/gIrO1MedhlsPsTSFljbO7aga+HM3NHNgvgVl/D5Ngglu3O4UhBJQOC2/hlAXy4JZ2Cihq+S8plV+Ya/nP9aEb379pnoUZzVrPsIcjZ1bl99hsBc9rOXtXIwYMHWbRoERMmTKCoqIiAgAAaGhqYMWMGiYmJjB49uim42po1axg+fDibN2+mvr6e8ePHd+6YTxG7bPoiYhWRHUAehuLeCMQA14rIFhFZJiKxZvUwoHls4gxT1pb8+Hvdafa5JT8/v8MTOoFRNxh/o8+D2Fnt1x94ATh7QuKHxj7AyahvZaN231J4qj/TgiuBk7tu1jfYeGd9KhMHBPLxXRMBuPqV9XyxPbP9cWo0mm4nMjKSCRMmAPDRRx8xZswYRo8ezZ49e0hKSsLJyYmYmBj27t3Lpk2b+N3vfsfq1atZs2YNU6ZM6eHRG9gVWlkp1QDEi4gf8LmIDAdcgWqlVIKIXIGR8Py0Z6WUWggsBCP2zun2R/AguOI16D8B7El87uwOo66DLa8bJ3Uv/oexKdxykLD4V7B/GdyzCTybbfCufxHqqwnN+YFw/6GsOVjATROjWr3V93vzyCw5yp/nDmV0f3++vm8Kd769hQc+TWRgHy8d90ejaQ07VuRdRWMI5SNHjvDss8+yefNm/P39ufXWW6muNs7rnHfeeSxbtgxnZ2cuuOACbr31VhoaGnjmmWd6bNzN6ZD3jlKqBFiJYZbJABpzfn0OjDTfZ2LY+hsJN2VtybuekVeDX0T79Rq56Bk4/xHDvXPhNMje2bJ81d9g+ztQVQCbFh6T5+2DVCOpgiR/x5TYINYfLqS+jZAMb/10hDA/96a9AF93Z166cQyBni788t1tlFZpG79GcyZSVlaGp6cnvr6+5ObmsmzZsqayKVOm8O9//5uJEycSHBxMYWEh+/fvZ/jw4T044mPY470TbK7wERF3YCawD/gCaDR6TwUOmO8XAzebXjwTgFKlVDbwLTBLRPzNDdxZpuzMw2KFqQ/ALV9BbaWh+L+4B0rSjU3hH5+G+Pkw+CIjpEONueWx9U2wOBsmpZR1nBfpTnl1PTvSS064xb6cMjYcLuKmiS2TtQR6ufLijWPILj3K7z/egU179mg0ZxyjRo1i9OjRxMXFccMNNzRl0wIYP348ubm5nHfeeQCMHDmSESNGtJnDo7uxx7wTAiwSESvGQ+IjpdQSEVkLvCsiv8XY6P2ZWX8pcBGQDFQBtwEopYpE5HFgs1nvscZN3TOWqMlw90+w5h+w6VXY9RHYGmDA+XDJvyFrh+EOum0RjL0Vdrxv+PrH3wA732OacxIuTi58vSubhONi+C/6KQU3ZwvXnXPiL5Ax/f3508VDWbB4D796fzt3TIlmdITfGfOPRqNxRI4PuPbWW2+1Ws/d3Z2ammP7fQsXLmy1Xk9hj/dOIjC6FXkJcHErcgXc00Zfb2DY/s8ePALgwidh/F3GCr8sC65+C6zOEHGOEePnpxeME701pXDOHRB+Drj64J6ygumDb2FJYjZ/unhoUyiIkqpaPt+eyWXxYfh5uLR625snRpJfXsOb647w9a5s4vp588DswUyP69uNk9doNL2N3nsit7Pxi4B5L8BNn4GbzzH5lN9CeRZ8+ycIHgL9JxoPhJjz4eB3XDIyhPzyGjYeORaL5531qVTX2bhtUrQhqCyE4w56iAh/uHAwGx+5gCcvH05dg407397K90m53TFbjUbTS9FK/3SJmQH9RkL9UWOV32iCiZ0F5VnMCMjHw8XKV2aI5uq6BhatT2Ha4GAG9/OG0kz493AjFHTd0RO693J14sbxkXxxzySGhvrwy/e26QieGofkTM7y11Ocymeilf7pIgIXLIDQMTDymmPygTMBcDvyHTOH9mXZ7mxq6218ti2Tgopa7jzPPJ23+TVD2SevgHevPrYpfBzebs68ffs4BgR58rNFW9iaemZvh2g0nYmbmxuFhYVa8TdDKUVhYSFubm4daterc+T2OP+dCk5urDj3be5YtIXXb0ngya/34unqxOJ7JyH11fDPoRB5rrEB/PldEDYWrn4TfMNb7TK/vIbLXlxHiK8bn9x9bjdPSKPpGerq6sjIyGjyhdcYuLm5ER4ejrOzcwv5yXLk2nU4S3OKxM6CNc8yJcyKj5sTCxbvIaP4KP+5frThiZP4ERwtMhKzR002DoZ9fBv8axgExED0FJj06xZhoIO9XZk2OJivdmahlNIePRqHwNnZmejo6J4eRq9Am3e6ksGzQdlwSV7GnOEhZBQfJdzfnTnD+xmneje+An1HQKTp4zvkEvjlBpj1JATFws4PYfF9J3Q7INiLsup6iipru3lCGo3mbEcr/a4kdIyxYk/8kHnxRkC2n08ZYBzGOrIa8pJgwl0tw0MEDYRz74UbPoQpvzcCvxWntOh2QJB5FLygsrtmotFoegla6XclIjDyWkhZy7nB1Xx810RumhBplG14GTyCYPhVbbcfdR0gxqGvZgwINpT+4Xyt9DUaTcfQSr+rGXkNoGDXx5wTFYDFIpCzGw58Awm3GTH828IvAgZMhR3vtfDjD/Nzx9kqHNYrfY1G00G00u9qAqIhYrxhn1fKeC39g3HSd8Iv228fPx9K05oCuYGR2CUy0JPD+a27d2o0Gk1baKXfHYy8BvL3Qu5uI05/2nq44C+G4m+PIXPB1Qe2v9tCHB3kqW36Go2mw2il3x0Mu8KIvrlpoZGcPSwB4m+0r62zOwy/ApK+hOqyJvGAYE9SC6t0fl2NRtMhtJ9+d+ARYPjsb3sbELjxE7B04HkbPx+2vgXf/tHoqySNGdbR/LdhIJnFR+kf6NFVI9doNL0MrfS7i5HXwP6vjfg8ofEdaxueAH2HG4lbrK4gFoZ7JQELOFRQoZW+RqOxG23e6S7i5sIlz8EF/9fxtiJw+zfwu33wSA4k3I57RTqgONIFbpsZxVU88+0+8sq758j7/pxynvv+oI6rotF0A3ql311YnYxEK6eKq/exXL0B0Uh9FdFulRwu6DwPnroGG2+uO8K/vjvI0boGKmsa+L9Lh51Wn3//Zh97ssr49QWxjOnv32qdF1cms3hnFpeNDiUy0POU76WU4mBeBVGBnrg46fWMRtMa7Sp9EXEDVmMkQncCPlFKLRCRtzDSJJaaVW9VSu0QIxjMcxjZs6pM+Tazr1uAP5n1n1BKLerMyTgM/lEAjPMr7TQPntKqOq5/dQNJ2WVcMKQPDTbFp1sz+MOFg/FyPbW1waYjRby06hDOVuHHA/nMGtqXh+bEMSDYq6lObb2NlfvyANiaWmyX0s8rq+am1zcR5u/O+OgARoT5suFIEV9szyStqIo7Jkfz57lDT2nMGk1vx57lUA0wXSk1CogHZpu5bwHuV0rFm68dpmwOEGu+7gReBhCRAGABMB4YBywwc+VqOoq/EXhqpEdxp53KXbwzk6TsMp67Lp5Xb07gVzNiKa+p5/Ptp5a7vrbexiOf7yLMz52fHprB72YO4qdDhVy7cAM19Q1N9TYcLqS8ph6ALanFdvX9+fZM9ueWk1JQyd+W7eOG1zbynx8OEhHgzrkxgbyzIZXs0hNzE2g0GjuUvjJotCE4m6+TGV/nAW+b7TYAfiISAlwIfKeUKlJKFQPfAbNPb/gOil8EIAx0yie7tJqq2vrT7nLprhwG9vFiXnwYIsLoCD+Gh/nwzvqUU7K1v7rmMAfzKnhs3jCCvV25b0YsL944hvzyGr7ZndNUb3lSDu7OVsZFB7A1xT6l/+WOLOIj/PjhD9PY9MgM3rg1gfUPzeDdn03g6StHopTihR+SOzxmjcYRsMvwKSJWEdkB5GEo7o1m0ZMikigi/xIRV1MWBqQ3a55hytqSH3+vO0Vki4hsyc/P79hsHAUnV/ANJwzDLHK6Jp6CCiOd40XD+zXJRISbJ0ZxILeCDYc7lrAlrbCK51ccZPawfswYciyn75SBQUQFevDO+lQAbDbFd0m5TB0UzOSBQRzIK6f0aN1J+z6YW05SdllTALs+3m5Mj+tLP18jnEVEgAfXnhPBR1vSSS+q6tC4NRpHwC6lr5RqUErFA+HAOBEZDjwMxAHnAAHAg50xIKXUQqVUglIqITg4uDO67J34RxFQY5heGpV+cWUtqYUdfwAs35OLTcHs4SEt5JeOCsXPw5l3NqSctL1SiieWJDHznz8y+ekfuOj5NThZhAWXDoWiI5D4Mez6BMveL7h3eD1bUotJyipjZ0YJuWU1zBrWl7GR/igF29NOvtpfvDMLi8DFI0ParHPv+bGICM+vOGj3Z6DROAodcnFQSpUAK4HZSqls04RTA7yJYacHyAQimjULN2VtyTWngn+U6bZpRNvccLiQmf/6kateWd9hc8yy3dlEBXowJMT7mPBoCW7OVq5NiODbPbkntZEv253Da2uPEOTlyrjoAObFh/LqzQmEeDnBu1fBZz+DT++Aj2/lisQ78XKy8b+NqSxPysVqEWarn5iw4mpcLDa2nsSur5Tiyx1ZXBwt9PngYlh6P2TtMOIZNaOfrxvzx0fy2fZMHZ9IozmOdpW+iASLiJ/53h2YCewz7fSY3jqXAbvNJouBm8VgAlCqlMoGvgVmiYi/uYE7y5RpTgX/KKQylwE+8NGWdG58bSNlR+vJL68hr7zG7m6KK2v56VAhc0aEGFm4MrfC+9fD05Gw4z3mT4jEphQfbk5vtX15dR1/+WoPw0J9eOeOcfzzmnievHwE5w4MMg6TFSbDvJfgnk1w6QtYjhZx/4BUvtieyZLELCYMCMBj83+wZm/jwuAitpzErr89vYS0oip+5fYNZG2DrYtg4VR4ZTKseAwO/QC1xi+du6fFYLUI721M69jnqtH0cuxZ6YcAK0UkEdiMYdNfArwrIruAXUAQ8IRZfylwGEgGXgV+CaCUKgIeN/vYDDxmyjSnQqPbpn85GcVHuWBIH166cQwA+3LK7e7mu725NNgUcwd5GYnZX50OqT+BXySseJwIb+HcmEA+3ZaBrZU4P/9YfoC88hqevHyEkRymkdpKWPUU9J8I8TdA8GAYdT149uEyWUNVbQPpRUe5LrwYchIBuMAnkx3pJdQ32E64D8DiHVkEOVURm/EJDL8S/rAfLv4HuHjBuufgncvh6SjYt5Rgb1eGhPiwJ6us1b40GkfFHu+dRKXUaKXUSKXUcKXUY6Z8ulJqhCmb3+jhY5p87lFKxZjlW5r19YZSaqD5erPrpuUABBhumz8fLjx79ShemT+WsZGGB+yBDij9ZbuyCfd3Z2jGB3BwOUz/M/x2N1z2MpRnwcb/cvXYCNKLjrIppeUzeldGKW+vT2H++EjiI/xadrz+JajINaKJNmYGszrBiKvxTV/BpDArAOdXf2eElnDxZpT1CEfrGtibfeL46xtsLEnM4s99NyC1lUbuYHd/OOdncMe38GAqzP8UfMKMBwAwNMSbvTll+qSvRtMMfWzxbMX01Y9xKuCqseGICP6eLvTxdrV7pV96tI61yQVcMiwQ2bQQYqbDeX8wTv5GTTKCxK39JxcOcMPb1YmPt2Q0ta1vsPHIF7sI9HLl/tmDW3ZcWWAo3ri50H98y7JR14Gtjr8NOsCC2QPw2vepET46bAyhlXsB2JLa8uFSW2/j0cV7KK+oYHbFFxAzA/qNaNmvqxcMvADG3gLpG6DoMENCfCipqiO7tHvCSWg0ZwNa6Z+tuPuDq6/hHdOMwf28OZBrn9L/dk8OdQ2K69w3Gavyife0rDDjUaguw33T88wdFcKy3dlUmgep3lh3hMSMUh6dOxQfN2eoOwppG2Hz68ambV0lzFhw4k37jYA+w+ifvpjbgvZCdQmMng9hY3Ep3Eukj6XFZm5eWTXXv7qB9zam8e8h+3CtKTBW+W0x4hpAYOeHDAnxAWBvtjbxaDSNaKV/tiIC/pEnJE0f3NdQ+vbE2f90awbRgR703/8G9BlqrKCb028EjLgaNr7C9XFOVNU28PWubA7nV/CP5QeYObQvc0eGGAr/v+fBG7Pg699B5nY4/48QPKj1cY+6DjI2w49/B98IiJ4KYWPAVs8lfQvZeKSI19ce4U9f7OLi/6wlKauMF64fyZzSjyEkHqLPa3tSvmFGeeIHxPU1wj1opa/RHEMr/bOZgGgobrnSH9TPm5p6G2ntHExKL6pi45EifjMgE8lLMlb5jbb35kx/BBrqGJH1MQOCPPl4SzoPfboLVycLT1w23PD4Wf0MFByAS56H3+yGh1LhvPvbvvmIq0EsRjax+BvAYoVQYxN6qlc6+eU1PL4kicU7shgQ5Mnnd49jbtqzUHQIJv+m9XE2Z9T1UJyCd/42+gd4tLpHoNE4KjrK5tmMfxTsXwa2BkNxZu9kfGUiEMD+nDKig9oOXvbptgxE4MKyT8Crr6GI27pHzPnIrk+4csz1PLPcOPD096tG0tfHDfL2wbrnDUU79hb7xu0TAgOmGS6W8TeYslDw6stYpyMs+dUvCfF1I8DTBak7Cp/cDgeWweTfwdDL2u9/yCXGL46d7zMk5CaS9Epfo2lCr/TPZvyjoaEWyrKgoQ4+vo3+a/6As9SzP6ftQ0k2m+LTbRlc278ct9SVMO7nRmiHthh5LZSmcV3fDCwCU2KDuHpsONhssOQ3xibqrCfabt8as540/PdN11NEIHQMluztDA/zJdDLFakuhUWXwIFv4KJn4YIF7a/ywRjPkEth9+cM7+tKSmFlp8Qn0mh6A1rpn800KsziFNi2CIoOIfXVTPErYn9u26vbzSlFpBcd5T750Ei6nnDHye8TdzE4exJ46As+vmsiL1w/xjDr7PifkeR95mPgGdSxsfcdCqOPyxMcNgYKDh7LBfzD45C1Ha59x3gwdYRR10JNKVNsm1GqY2cXNJrejFb6ZzOmrz65u2HV0xAwAIDzvDLZfxIl98nWDCa5HiY05weYdJ+Rd/dkuHgaJpM9XzA21ANfD2cozYTlfzYPX83vnPmEjgEUZO+A3CTY8oaRXnLIJR3vK3oquPkxsNw4JqI3czUaA630z2Z8wsHiBD8+DZV5cPl/wcWbeOsRUgqrqK5rOKFJVW09S3dl8bjnx+DZB8bfbd+9Rl4DNaXGAS6bDT7/hWFauvSFjiV5Pxmho42/mVuNJPCuPjDt4VPry2KFPkPxLD+Mt5uTVvoajYlW+mczVifD5fFosWHDjhgHIaOIrDlAg01xqJVgY9/vzeOc+m0MqNoJUx8w7N/2ED3V2PBN/BDW/wdS1sCcpyFoYOfNxzPQCP+w6VU4vNJQ+O39CjkZwYORvL0M6etNkg7HoNEAWumf/fhHgViNg1QAofH4lu3HSkOrJp71B/N4yOVDlH8UjLHT2waMB8zwq+DAt7DiccPkMvqmTplCC8LGQlkmBMYapp3TITgOqktICK5nX055q7GDNBpHQyv9s50pv4fLX4GgWOM6JB5LQw1DrVnsb+VkriV5OXGkIuc/Ak4uHbvXyGvAVgeewYZPvj2eNB0l/Bzj74V/Bavz6fUVbISHOMczj6rahnbPLmg0joD20z/biZ7S8jo0HoDzfbNIPG6ln1dWTWjFLhqcnbAOndfxe4WMgqkPGTF5TsfscjLG3mJ49gyYdvp9BccBMMiaCQxkb3YZUSc5u6DROAJ6pd/bCIgBF2/OcUk9wbyzKaWIOEmn1i/m5H75bSEC5z8M4WM7abCt4OLZOQofwLsfuPrStyYVi2gPHo0GtNLvfVgsEDKSOHWY7NLqFopu4+EihlrScA0b2YMD7EZEIHgwToUHGBDsRZIOx6DRaKXfKwmJJ6jyIG5WW4twyHsOpxIihVj6De/BwXUzwYMhf19TIDqNxtGxJ12im4hsEpGdIrJHRP5yXPnzIlLR7NpVRD4UkWQR2SgiUc3KHjbl+0Xkwk6dieYYofFI/VFuHFDNFzsyqa23UVRZi1O+Ea+evo6k9OOgMp9RgcZGrg7HoHF07Fnp1wDTlVKjgHhgtpn7FhFJAPyPq38HUKyUGgj8C3jarDsUuA4YBswGXhIRa2dMQnMc5iGnK0MKKKqsZeX+PDYdKSLOYuaL7TusBwfXzZibuaNccwE4kKsTpWscG3vSJarGVIiAs/lSpsJ+BnjguCbzgEXm+0+AGWby9HnAB0qpGqXUEYwcuuM6YQ6a4zE3c+NsyQR7u/Lxlgw2HSlimDUd5R5gbHA6CqbbZowYZq6OpJLUaHojdtn0RcQqIjuAPIzE6BuBe4HFSqns46qHAekASql6oBQIbC43yTBlx9/rThHZIiJb8vPzOzgdDdC0mWvJ2ckVo8NYuT+P7/fmMsY1C+k7rGv8689UfMPBxYuAqsO4OVtaPbug0TgSdil9pVSDUioeCAfGich5wNXAfzp7QEqphUqpBKVUQnBwcGd37ziEjoacXVwzyp8GmyK9qIKohhTHsueD8YALGoSlYD+D+nqfNBCdRuMIdMh7RylVAqwEzgcGAskikgJ4iEiyWS0TiAAQESfAFyhsLjcJN2WariBuLtRXE5O/kvgIP/pLHs62asey5zcSHAf5ptLXK32Ng2OP906wiPiZ792BmcBWpVQ/pVSUUioKqDI3bgEWA41BXa4CflBKKVN+nendEw3EAps6dTaaY/SfYCRZ2fEud543gOl+eYbcIZX+ICjPZkQQ5JfXUFRZ29Mj0mh6DHtW+iHAShFJBDZj2PSXnKT+60CgufL/HfAQgFJqD/ARkAR8A9yjlDox9q+mcxAxUhGmrOGi8FoWjFNGXlrTm8WhMOc80vTg0SYejSPTbuwdpVQiMLqdOl7N3ldj2Ptbq/ck8GQHx6g5VUZdByufNMIh5+4xvXo8enpU3U+jBw8ZQF/255QxMSawZ8ek0fQQOuBab8avP0SfBzveBaWMdISOiF8kOLnhXZaMn0c4+7WvvsaB0WEYejvxNxo5dEtSHdOeD0YWrcBYpPAgg/t6sz9HB17TOC5a6fd2hlwCLqb1zdHcNZsTFAsFBxjcz5sDuRUYvgUajeOhlX5vx8UThl5mvO8ztEeH0qMExUJJGnFBLlTU1JNZcrSnR6TR9Ajapu8InP+wkVzFr39Pj6TnCBoEysZIj0IADuSWE+7vgJvaGodHr/QdAd9wGPdzxwq/cDyBxjGSaDGihuzTbpsaB0UrfY1jYCp9z7JDhPq6sSfz5Ju5DTZFg06krumFaPOOxjFw9QKfcChIZsaQOXy4OZ3csmr6+rg1VSk9WseyXdmsOVjAukMF9PNx45vfnNeDg9ZoOh+90tc4DkEDoeAAP58ygHqbjTfWHmkqarApbn5jEw99tostqUUEe7myL6ecmnp9aFzTu9BKX+M4BA2CgoP0D3Bn7shQ/rchldKqOgD+tyGVnekl/P3KkWx4eAY/P28AALmlNT05Yo2m09FKX+M4BMZCbTlU5HLX1Bgqaxv438ZUckqreebb/UyJDeLqhHBEhDA/dwDt2qnpdWibvsZxCIo1/hYcYGj0eUwbHMwba4+wLbWYugYbT1w2HDE9nEJNpZ+llb6ml6FX+hrHIWiQ8bfgIAB3T42hsLKWFfvyuG9GLJGBnk1VQ3yNDV6t9DW9Db3S1zgOPqHg7Nmk9MdFB3BuTCClR+v4+ZQBLaq6OVsJ9HQhq7S6J0aq0XQZWulrHAeRJg8e41J487ZzEAQXpxN/9Ib6ueuVvqbXoc07GsciMBYKDzZdujpZW1X4AKF+blrpa3od9qRLdBORTSKyU0T2iMhfTPnrpixRRD4RES9T7ioiH4pIsohsFJGoZn09bMr3i8iFXTYrjaYtggZBSTrUnUSZKwXLHuLBvAfJKjmqI3JqehX2rPRrgOlKqVFAPDBbRCYAv1VKjVJKjQTSgHvN+ncAxWbO3H8BTwOIyFDgOmAYMBt4SUSsnTkZjaZdggYCCgoPtV3nx6dh48tEl2+lvvYoZUfru214ncmRgko+3pKuH1qaFrSr9JVBY6ohZ/OllFJlAGL4uLkDjf+y5gGLzPefADPMOvOAD5RSNUqpI0AyMK7TZqLR2EOTB8+B1su3vAmr/gb+0QiKCMk763z1bTbFW+uOMOe51dz/SSI7M0p7ekiaMwi7bPoiYhWRHUAeRmL0jab8TSAHiAP+Y1YPA9IBlFL1QCkQ2FxukmHKjr/XnSKyRUS25Ofnn8qcNJq2CYgBBAqTTyw78C18/TsYOBMu/y8A0ZJzVtn1CypqmP/6Rv7vqyTOiQpABH7cr/8faY5hl9JXSjUopeKBcGCciAw35bcBocBe4NrOGJBSaqFSKkEplRAcHNwZXWo0x3DxAN8IyNt7Ytmqp4xonNcsgmDjF0Gk5JJV2nlKf1dGKZuOFJGcV05RZW2n9dvISysPsTmliKeuGMHbt49jZLgfqw7kdfp9NGcvHXLZVEqViMhKDJv8blPWICIfAA8AbwKZQASQISJOgC9Q2EzeSLgp02i6l8hz4cA30FAHVmdDVpYNWdtgxqNGtjE8Ue7+xDTkklrSOb76uzNLufTFtTQ3sS+4ZCi3TYrulP6VUny3N4fJA4O4bpyRMGfqoGBe+OEgxZW1+Hu6dMp9NGc39njvBIuIn/neHZgJ7BeRgaZMgEuBfWaTxcAt5vurgB+UsZO0GLjO9O6JBmKBTZ04F43GPobOg+oSOPLjMdn+pcbfwRc3iSQghljnvE4z7/z7+wN4uzrx5q3n8Nx18SRE+vPCD8lU13VOJM+DeRVkFVVwTWg+HPoB9nzBJZ5J2BSsSS7olHtozn7sWemHAItMTxsL8BHwNbBGRHwAAXYCd5v1XwfeEZFkoAjDYwel1B4R+QhIAuqBe5RSOm6tpvuJmQ4u3pD0JQy8wJDtXwoBAyB48LF6AQOIzFnVKUp/R3oJ3+/N4w+zBnF+XB8Agr1dueHVjXy+PZPrx3UsleXLqw6xLrmAd+4Y1xQv6Pu9udxr/YI56z+F9Ua9WCDe/d/8uD+fS0eFnvY8NGc/9njvJCqlRiulRiqlhiulHlNK2ZRSk5RSI0zZjY3ePEqpaqXU1UqpgUqpcUqpw836elIpFaOUGqyUWtaVE9No2sTZDQbPhr1LDBNPdRkcWQ2DL2qZUjIwhsCGfAqKT9/75Z/fHcDfw5lbm5lyJg4IZFioD6+tOYytg1m6liRmsTa5gFXNNmlXJWVxs8tKiJoCt30DN30BwI3Bh/nxQH6H76HpnegTuRrHZOhlcLQIUtZC8vfQUAtxF7esEzAACwrninTqG2ynfKstKUWsPpDPL6bG4OV67Me1FKdw56RwDuVXdmiztaq2vinH739XG+cNCitq8Mv8gUBVBBN+CZETIeZ88IvkXMseCipqSMo+eYpIjWOglb7GMRk4A1y8IOkLw7TjEQgR41vWCTCCsEWSQ275qSVTqaqt59nl+wnycuHmiZHHCrYugv+M4ZKUvxLi68arq4+03clxJGaUIrY6pke7s+FwEYkZJfywL48bLD9Q59kPYmcdqzxgKiHFm7HSwI8HtOumRit9jaPi7A6DLoS9X8HB5TBoNliOOyDeqPQ76KtvsyleXX2Ya/67nlF/Wc6Gw0Xce/5APFycjBAPq5+Br+4DjyAsuz7ioWHFrD9cyO5M+8xI29KK+ZPT/3i19E7iXIv47+rD7EhM5DxrIk5jbwZrs6266KlYasqY1yePVfu166ZGK32NIzN0HlQVQnWpYc8/Ho8AGlx9iZLcDin97/fm8uTSvZRX13P7pGjeuWMct5wbBTUVsPR++OEJGHEN3LsZfMKZm/FvfFwtvL7WvtX+ttQSJrscxFqVzzse/2DtrkNEpHwECDL2lpaVo6cCcJlfMtvSSpj/2kYue3Edl7+0jtwyHTbaEdFKX+O4DJwJzh7g5G549LRGQAxRkkNWB3z1V+7Px9PFypf3TOLhi4YwJdIDWfccPDcSNr8KE+4xTvy6+8Gsx7Hm7eKv/bewJDGLvPKT30cpxZ7UXKJtaRA1haCadF50fo4rZBVFoVPBN7xlA69g6DucBNsuhoX6UFlbj4vVwva0EtYfKrR7Tpreg46nr3FcXDxgwt3GJq6LR6tVrEExDMhexXI7V/pKKX7cn8fN4Tm4rHsWMrdC2gbjXEDMdJj2MEQ0Czk17HLY8gZzsl/jTw2xvL8xnV9fENtm/2lFVQQdPYzVtQHO+Rky8homL/4VALWT72y9UfRUPDa/xuKHxoKzO7XVVTz4xBMcyo2ilUgoml6OVvoax2bGoycvDxhAPz4lz063zYN5FUSUbefBmschS4wAb3EXw5hboP/4ExuIwJy/Y31lMn/ts4L/2xjE3dNi2ozxvy2tmOEW0wwUMgoCoqkrzcZ2eDWucW1EKx8wFTa8COkbIXoqLkt/y7+cPuA/KeEYQW81joQ272g0JyMgBis2GopS7Kq+an8el1nXYnP2hAcOw72b4LKXWlf4jfQdCrGzmF6/lrzyapbtzm6z6va0EuKdUlFuvuAfBYDz+Q/iesfXJ25ENxJ5Llic4PCPsOElSPwAAGtxK0HnWsFmU9z3/naWJGbZVV9zZqOVvkZzMkwPHtfyVLuqr9mXzcXOW7DEXQQeAfbfZ8gluFVlcaF/Nm/9lNJmtW1pxSS4pCEho1oeJDsZrt4QNhZ2vAfL/wRxc6m2euFTlUadHecPfjyQz+KdWfxnRfJJY/P/bdleHv1yt31j0vQYWulrNCcjMAaAfnWZ7cbIqaipxyVtNT6qHIZd0bH7DJ4DYuWevklsTythZ3rJCVWqaus5mF1MZH2KYdrpCNFToSIHggbD5a9Q5RVJf3JILaxqt+lra41D9ftzy9mT1foBr8KKGt5Ye4R3NqSSUdx+n5qeQyt9jeZkuPtT6+RNpORSXHXyUMjrkguYI+upd/Y2Dn91BI8AiJ7C8NJVeLpYeGxJEqVVdS2qJGaUEq0ycVK1EBLfsf5HXGWEZ7j+PXD1RgJjiJYckvMqTtpsT1Yp65ILm/YZPtma0Wq9T7dlUNdg/Ap4f1Nax8am6Va00tdoToYIR70jiZKcduPfr9mbySzLFmTIJeDk2vF7DbkUS/FhXpzpTmJGCVe+8hPpRcdWzSds4naE4MFw65Imc5VnyGBCpYAjuUUnbfb62iN4uFi567wYZg7tyxc7Mqmpb/mLRynF+5vSOSfKnxlxffhwczq19acetkLTtWilr9G0Q73fgHaVvlKK2n3L8ZEqrCOuPLUbxc0FhGkNG3n79vHklVVz+Us/8bdle7nipXX8Y/kBJnlkGOEjAmJO7R4mLn1isYqiJLNl2sjmQdlyy6r5amcW14wNx3fpL/i190pKqupYua/lyd71hws5UlDJDeP7c+OESAoqavlmT85pjU/TdWilr9G0gwREESYFFFW07at/MK+Cc2tWU+3sZ7hIngrefaH/BNi7mIkxgXx697m4OVt4bc0RFHDX1AFcGJAL/UaA5TT/65p7FfX5xzx4CipqGPXYci56bg0LVx/i+RUHabApfjG4CnZ/yqCtj3G951Y+2doy99F7G9PwdXdmzvAQpsYGExHgzv822Lfxrel+tJ++RtMOrr59sIqisqQQaD3u/dqkdK6xbKNh8FXHsnGdCkMuhW8fhsJDxPaNYdVvJlKLEx6uzmBrgL8lwZibT73/Rhq9kkqPYLMpLBbhm905lFfXo4C/LjVyIs0Z3o+Q7G9BLBAyiseyX+DaA37kl48g2NuVgooavt2Tw00TonBzNlxGbxwfyVPL9rE/p5zB/bxPf6yaTkWv9DWadnD3NZKe1JS1HaWyZPc3eEk1nqOvOr2bDZlr/P3kdnhhHE5PheLx+lQoSTOSuddVdXwTtzU8Aqhx9iXMlt2UA/j7Xaks8nyBpVd58sPvp/LIRUN45OIhRt6BiAlw4yco73781+kZnv1wOW+sPcI/lh+grkFxw/hjmVCvHhuOi9XCuxv1av9MxJ50iW4isklEdorIHhH5iyl/V0T2i8huEXlDRJxNuYjI8yKSLCKJIjKmWV+3iMhB83VLW/fUaM4kLJ6BANSWt55ysLKmnn5566ixeEDkpNO7mV9/I2VjVaGxGp94L5RmwKszjHDM0PFN3Dao9Y0myvTgKaqsRaX8xNSGn5AVf2FAsBc/P28A4SoH8vYYDyPPIFxu+gRPq41b0x/h70u28/6mNMZHBzCwz7EVfaCXKzOH9WXpruyT+vVregZ7zDs1wHSlVIWp2NeKyDLgXWC+Wec94GfAy8AcjCxtscB4UzZeRAKABUACoICtIrJYKVXcmRPSaDodD0Pp2ypbV/obDhUwWXZSEXourk6dkHz8+vdaXo+eD+9eZYRScHIzQjt0As59BhKVv5JleRXkllUzUXYZBYdXQdYOCI03VvlgbjIDwYPwuO5N4t67hp1jv2H/hKeJCPA8oe8J0QF8nZhNRvFRIgJaj2uk6RnsSZeolFKNzrzO5ksppZaaZQojwXljeL95wNtm0QbAT0RCgAuB75RSRaai/w6Y3dkT0mg6HXfjZK0cbd29MXHXdvpb8vEZPqvV8tMmeDD8bAX0PxdiZ7aMl38auPaJJUwKSc0tZOmuHKY7J6FCRoGrD6x7zqi0b4mxcezfLAHMoFnI1Adx3fMhI3M/x9/zxAfd6P7+AGxv5ZCZpmexy6YvIlYR2QHkYSjujc3KnIGbgG9MURiQ3qx5hilrS378ve4UkS0isiU/X2f60ZwBmCt9a/WJP0qVUqjkFQA4x17QdWPw6gO3L4Nr3um0LiVwIAC5KfvYe+gwg9RhJO4SSLjNyCiWthHSN0HcJSc2nvqgEZp66QOw6xMoOAh1x7ybBvfzxs3ZwvY0/UP+TMMupa+UalBKxWOs5seJyPBmxS8Bq5VSazpjQEqphUqpBKVUQnBwcGd0qdGcHi6e1IkLLjUlJxSlFFYxonor5e7hTW6QXYq98XbswfTgkaJDJKg9hmzANBh/txGg7eNbAHVsc7k5FgtcsRB8w+DTO+CFBHiyH3zxSwCcrRZGhvmxPa2k88ar6RQ65L2jlCoBVmKaZURkARAM/K5ZtUwgotl1uClrS67RnNmIUO3ki3t9SYvDSwCr92Yy0ZLU8bALZwLmQypKcpjlthfl6gOho8EnBEZeC+XZ4B8NfYa23t4jAO5aB7ctg8sXGvkCdn0CtcYp4tH9/UjKKjvhBK/GDiryoKGu/XqngD3eO8Ei4me+dwdmAvtE5GcYdvrrlVLNz1wvBm42vXgmAKVKqWzgW2CWiPiLiD8wy5RpNGc8ta7++FFOeXV9C3n27h/xkmq8h7URy/5Mxs2XWtcAoiSH85ySkKgpx/YLzr0PEBh66cl/Xbh6GaGbR11reBo11ECK8aN/dH8/ahtsbQZp05yEr34Nr7aRze00sWelHwKsFJFEYDOGTX8J8ArQF1gvIjtEpDEbxVLgMJAMvAr8EkApVQQ8bvaxGXjMlGk0ZzwNbv74SzmFlTVNsuq6BgKy19CA1QhmdhZiCRrIBa5JBNRmGqadRoIHwc9/gPPut7+zyElG+smD3wHNNnO1iadjVJdB8oou+zfVrhuAUioRGN2KvNW2pjfPPW2UvQG80cExajQ9j0cAAaS1iLS5OaWIc9lJefBo/Nx8enBwp45TcCzBmZuMi+ZKHyBszAn1T4qzm6Gokr8Dpejr40aorxs7tAdPxzi43PjFNHRel3SvT+RqNHZg9QzCTyooqjxmZ01NTWG4JQXXuJk9OLLTxNzMxTsUgtrOzWs3sTOhOAUKDwHGal978HSQpC/AOwTCz+mS7rXS12jswNk7CD8qKK44FurYKX0dAO5ns9Jv9DgaMK1zPIMGmm6ryY0mHj8yio+SV159+n07AjUVhnlsyKWnH1SvDbTS12jswM0nGKsoKkqPbUNZSlKMN8FxPTOozqCPmRi9s84YBERD4MBmdn0/AHZou759JH8H9dVdZtoBrfQ1Grtw8THOjDQPuuZSmU2FxQdcTgxDcNYQPAjuXt/x9I4nI3YWpKyF2iqGhfribBV9Mtdekr4Ezz5GiO0uQit9jcYezCTn9WbQtfoGGz61eVS69e3JUXUOfYd27qGvgReYrptrcXO2MjTER9v17aG2Cg4sNw7DWaxddhut9DUaezDj7zQGXcsurSaEAuq9Togkoml03TTt+glRAWxLLSGvTNv1T8qhFVBXCUPnkVJQSXE76TlPFa30NRp7MOPviBl/J7WwihApxOoXfrJWjomzG0RNNqJ1AjdNiKTeZuO1tUd6dlxdRW0VVBZAWRaUn0aayN2fGYuLyMk8ungP17+6ofPG2Ayt9DUaezCVvlO1sZGbmZePn1TiERx5slaOS+hoM+lLNVFBnlwyKpT/bUjtstVrj5G7B56KgGdi4J9D4B+DjYNVHSVtA+z5HOJvwCZWtqcVMybSv/PHi1b6Go19uHhSL8641pYAUJpjrFq9+kT13JjOZPoMAWWDAiPx+j3nD6SqtoE31/Wy1f7B5WCrhwv/CnP/Zcgyt3Wsj7pq+PJe8IuAaQ9zKL+C8up6xvTXSl+j6TlEqHb2w7PBCCB2tCANAIs277ROY5C2vL0ADOrrzYXD+vLmTymUVXdNILEeIWUtBA2GifdAwu3g2x/y93Wsjx+fgsKDcMlz4OrF1lTDhDjGdHftbLTS12jspM7VnwApp7iyDlWaYQh99EZuqwQMAKuLkWrR5N7zYymvrued9b0kd25DvWGWiWqWIjN4MBTst7+PrB2w7nmIn29EKQW2pRXj7+FMdFDXuAJrpa/R2InNLaAp6JprZTYKAZ/Qnh7WmYnV2UjraK70AUaE+zJ1UDBvrjvSO3Ln5uyE2gpeTg3hT1/s4qMt6RR7RhsJZWx2hJM++D18dBN4BsGFTzSJt6WVMKa/P9KZbrTN0Epfo7EXjwACKOdQfiVBDXlUuQYZyk3TOn2GtFD6AOfGBFJQUUtFTX0bjc4iUowwHG+kh/Hl9iwe+CSRp7Yo40RtSVrb7QoPwbvXwLtXgsXZyIbmbtjvS6pqSc6r6LJNXNBKX6OxGycvI+ja9rRiQqRQ++i3R58hUJpuhAo28fcw8umWVPUCu37qOjIsYYT3j2LnglksvW8KB23mvwlzA/sECpJh4TRI/QlmPg6/3AD9xzcVN55c7qpNXNBKX6Oxm8aga4lphYRIkd7EbY/GzdxmG5uNSdSbh6g+K7E10JCyjtW1g7hkZCgWizA01IdC9yijPL8Vu35NBXx4o/Hr8O51MOk+cGqZVH5bajFWizAqwrfLhq6VvkZjJ26+RtC19KwcwqQA96D+PT2kM5smD56kJpG/h2EOKz7bV/o5u7DWlrNRDeXikSFNYv/APpRY/E9U+krB4l8ZvwCuegP8Wz/fsS2tmCEh3ni4tJvq5JSxJ12im4hsEpGdIrJHRP5iyu8VkWQRUSIS1Ky+iMjzZlmiiIxpVnaLiBw0X7d0zZQ0mq7B4mn8Mw+3ZeImdTj5a6V/UnwjwMWrhV3fr8m8c3av9FXKWgBqwyfS18etSR4Z6MEhwk704NnwMuz5DGY8emKyGpMGm2KHuYnbldjzOKkBpiulKkTEGVgrIsuAdcASYNVx9ecAseZrPPAyMF5EAoAFQAKggK0islgppSMxac4OzKBrwywpxrV21zw5FosRdrq1lf5ZfjK3fP+PFNn6MnnMyBbyyAAPkupCGJO/EVHKCGRXWQjfPQpxc2HSb9rsc39OOZW1DV2u9Ntd6SuDCvPS2XwppdR2pVRKK03mAW+b7TYAfiISgpFE/TulVJGp6L8DZnfKLDSa7sAMujZMUoxrX23Tb5fjPHh83XuBecdmwzljAxvVUOYMD2lR1D/Qk4O2UKSm7Fgcnn1fga0Opj100mim28xIpGO70HMH7LTpi4hVRHYAeRiKe+NJqocB6c2uM0xZW/Lj73WniGwRkS35+fnHF2s0PYcZf6dppa+Vfvv0GQqV+VBh/F92slrwcXM6q807KnMr7g1llPQ5hwDPlhuxkYEeJKtGDx7TxLPncyOxTN/hJ+13c0oRQV6uhPu7d8Wwm7BL6SulGpRS8UA4ME5ETj7600AptVAplaCUSggODu6q22g0HcdU+oMkgwaLM3gEtdNAQ58hxt/mJh5Pl7NjpX9oJbx/vXFqtpGs7dS/ey3FyovQsXNPaBIZ4HHMbTP/gBF988hqGHpZm6v8kqpafvfhDr7ckcW0wcFddiirkQ5tESulSkRkJYZZZncb1TKBiGbX4aYsE5h2nHxVR+6v0fQojUHXqKPGMxJrF+Uw7VU0j8EzYCoohb+705nvsllTDl/8Esqz4MA3MOGX0H8i6vNfUNzgwd3WJ1k0dtgJzYK9XalwDqTa6oVb/j7Yu9gIPDfs8hb1lFKkFx1l3aEC/rH8ACVVtfxq+kDunT6wy6fWrtIXkWCgzlT47sBM4OmTNFkM3CsiH2Bs5JYqpbJF5FvgryLSaLCaBTx8esPXaLoREWpd/HGqycPqr007duHVx9gLyUuCjC2w7EFeL07lNtsbPT2yY5SkQXEKRJ93TLbqKUPh3/Ax7FsC61+A9S9Q7hvH3Nxfcd9lU/ByPVF9igj9AzzJrI4gpuCAEUgtMBb6HntA/HXpXj7blklBRQ0Aw8N8ePv2cQwN9enqmQL2rfRDgEUiYsUwB32klFoiIvcBDwD9gEQRWaqU+hmwFLgISAaqgNsAlFJFIvI4sNns9zGlVBEazVmEm28w5OVpd017ETFW+7s+gW2LACEQRV1VSU+P7Bif3wWp6+C8++H8R4wY+RtehjG3wKBZxmvUddTvW8a1O8bj19eb686JaLO7/oEeHEgPIyZ7E9RWwJQ/NJl2dqSXsHD1YaYNDmbm0L6MjfRnUB9vLJauNek0p12lr5RKBEa3In8eeL4VuQLuaaOvN4Az6BGv0XQMi+m2qd01O0DEOMjYBJN/ZwSoW/oHqDr99d7W1CKGhPic3kGm3D2Gwg8aBKufgaIjxsrf3Q8u+L9j9SLP5Y3Uvuwt3sfbtw/Fydq2aS8ywIPEg32ZYzXDTzQz7bz9Uwperk68cMOYVn8pdAfaKKnRdARzMxdfrfTtZtrDcP8huGAB+BknUd3qSqitt51yl4kZJVz58npeWnno9Ma26VVwcoPbvzWU/O5PjAfUrCeazmUAFFTU8J8VyZw/OJjzBp3cwSQy0IN9DWb01aDBTZvZhRU1LEnM5ooxYT2m8EErfY2mYzQp/bZ/3muOw8kF3Ex7tfn5+Uv5abltPr8imVjJ4JsdKacepvloCSR+CMOvMhT85N/Cte8a5phR1zdVyyuvZv5rG6mpt/HIxUPa7bZ/oCcHbOaez7DLm0w7H2xOp7bBxs0TezbFplb6Gk1H0Oad08P8/AKk/JTdNpOyyijct5ZvXR9katlidmWW2tVu7cEC1iUXHBPsfB/qqmDcz47JhsyFGX9uUtRphVVc9fJ60oqqeOPWcxjYx7vd+0QGeJBJMD+OXwiTfg1AfYONdzekMmlgoF19dCVa6Ws0HSFoMLj6gp/eyD0lGlf6lJ+y2+aLP+zjCZdFWFDEWw+zJDG73TYpBZXcsWgzt7252UhHaLPB5tdQYefwyEYnzn92Fc99f5CskqMApBdV8cGmNK585SfKqut47+cTmBxr37mMMH93rBZhsyUeXDwAWLEvj6zSam6eGHVKc+5Mes6wpNGcjQy/EuIuApeuSWXX63H1xmZxxl8qTsm8cyC3HO+9HzLM+TC4+zOmPpOnErN5aHZcmx4wSin++FkiC61/Z7gcYe1bYxg8MQGvwmQ+j3yUdzemEdfPm399f4B/rzhAH29XcssMd8rIQA9eu3k8sX3tX507Wy2E+bmTWlTVJHt7fQphfu7MiOvT4Tl3NlrpazQdwWLRCv90EEG5+eNfe2rmnTe+28aDTh9SFz4B5wFTCF39DwrKS9meXszYyIBW23y0JR3nlJVMddnG0b5jmJ7zE14/raDSyZ+H9sfws8nRPHLxEDKKj/Lx1gyOFFRyTpQ/58YEEhPsdUonZCMDPUgrrARg1f481iUX8tCcuJN6/XQXWulrNJpuRTwDCSgv51AHV/p55dUM2/8ffK2VWOb+A4oOI9gY5pTJVzuzW1X6eWXVPPn1HhZ7fIzyjsL959+ybE8273z4ASU1Xlw5LoZHLh6CiBAR4MHvZg7qlDn2D/Dg613ZVNbU88jnu4kJ9uS2SVGd0vfpopW+RqPpVsQzkEBLPls7uNJPSznMDZbvyRk0n9B+w5vs5ZeFFPLCrmz+PHco1uNMPI8tSWJmw1qiOAzTXwcnF+aMiiSn4hqyS6t5cHZcl8S6iQz0oKSqjke/3ENW6VE+/sVEXJ2snX6fU0ErfY1G062IRyBBcqTDMfWPpm3DKgrLiCsNgV8UuHgxxTuXR9Nr2JxSxIQBgU31c0qr+X5XGht8PgX/kTDsiqay2yZFd8ZU2qR/gGEC/HRbBjdPjCQhqnXTU0/Q8wYmjUbjWHgE4icVHbbpq1wjUmfggHhDYLFA3+H0rzuMh4uVDzent6j/6bYMrrOswK8m2zgY1o0B8iIDjV8hob5uPDA7rtvuaw9a6Ws0mu7FIxAfVU5pZXWHmrmXHCCPQJw9myUZ6Tcca94erk0I56udWWSaLpdKKT7bnMJ9rl9D1BSImdGZM2iXAcGeTB0UzLPXjOrR07etoZW+RqPpXtwDsGCjrqpjmVIDqw6R5RrVUthvBNSU8YtRhmJ9bc1hALakFhNZsp4AWyGMv+ukGau6AlcnK4tuH8e5MWdezgWt9DUaTfdiHtCSjgRdszUQVp9Oqddx8eb7jgCgX9VBLh0Vygeb0imurOWjzenc4PwjyiMYBl3YWSPvFWilr9FouhdT6TvVFNsdN6c6NxlX6qgNPM4+3mcIiAVydvOLqTEcrWvglR8PsWHXXs6XbUj8dWB17uwZnNVopa/RaLoXM/6OjyqjrLreriZFKTsAcO53XLYqFw8j/2zubgb382Z6XB/+u/owcxpWYaUBRt/cmSPvFWilr9FouhdzpR/QgUib1Zm7sSnBt38r6bn7DoecRADunhYDKOa7rkZFjIfgzjls1ZtoV+mLiJuIbBKRnSKyR0T+YsqjRWSjiCSLyIci4mLKXc3rZLM8qllfD5vy/SKiDW0ajSPSIuiafW6blvx9pKtgwvoGnljYb4SR+ORoCedEBbBgVDn9bZnI6Js6c9S9BntW+jXAdKXUKCAemC0iEzDy5P5LKTUQKAbuMOvfARSb8n+Z9RCRocB1wDCMxOovmSkYNRqNI+Hiic3iQoBU2B1p07P0AMkSQbCX64mF/YzNXHL3QEUet1mXgYvXCcnINQbtKn1lUGFeOpsvBUwHPjHli4DLzPfzzGvM8hlinHOeB3yglKpRSh3ByKE7rjMmodFoziJEsLkH4I+d5p36GgKq08h1HdB6yIRGpf/etfBsLCR9CaPng6tX5467l2DXqQFzRb4VGAi8CBwCSpRSjbswGUBjVokwIB1AKVUvIqVAoCnf0Kzb5m2a3+tO4E6A/v11zHKNplfiEYh/WTkZlXaYdwqTsWKj3Ce29XKvvhA/H+oqIWwshI6B/hM6d7y9CLuUvlKqAYgXET/gc6DLzhUrpRYCCwESEhJOMQ+aRqM5k7F6BuIvOey2Y6WvcpMQoCG4jVSFInDZi507wF5Mh7x3lFIlwEpgIuAnIo0PjXAg03yfCUQAmOW+QGFzeSttNBqNAyGegQRZ7Iu/U521mzplxa3f4G4YWe/HHu+dYHOFj4i4AzOBvRjK/yqz2i3Al+b7xeY1ZvkPyjiBsRi4zvTuiQZigU2dNA+NRnM24RFIgJ0pE2uz9pCi+hEW5NsNA+v92GPeCQEWmXZ9C/CRUmqJiCQBH4jIE8B24HWz/uvAOyKSDBRheOyglNojIh8BSUA9cI9pNtJoNI6GewDeVNgVdM2pcB/7VTgxAR7dMLDeT7tKXymVCIxuRX6YVrxvlFLVwNVt9PUk8GTHh6nRaHoVHoFYUNRVthN0rbYSz8p0DtjGM00r/U5Bn8jVaDTdj71B1w58A0C6S/QZF6L4bEV/ihqNpvsx4+9YqgvbrlOWBV//nsPOsaT5TOqmgfV+9Epfo9F0P+ZK36O+lJr6Vrb2bDb4/BdQX8OjTr+hX6DexO0stNLXaDTdj7nS95dy0ouOnli+/j9wZDW22U+xsSyA/tqe32lopa/RaLqfxkiblLMvp6xlWeEhWPE4DLmEzKirqGtQRPhrpd9ZaKWv0Wi6H2cPlJMbAZYK9mYfp/RT1oCtjqxz/sjP3t6KRSA+wq9Hhtkb0Ru5Go2m+xFBPALpL0f5OLu8ZVneXhqcPLj4nTRsWHj79vEMDfXpmXH2QrTS12g0PYNHAGG2oyes9CvSE0muDSHY351Xb04gMtCzhwbYO9HmHY1G0zN4BBJsqSCrtJrSZjF4JH8f+1UE7/98glb4XYBW+hqNpmdwD8BHGav8vY2buZUFeNYVUeYVQ2BrCVM0p41W+hqNpmfwCMStrhSAfaaJpzZ7DwCe4SN6bFi9Ha30NRpNz+ARiFSXEOxhZa+5mZt5YBsA4XFje3JkvRqt9DUaTc/gE4KgmBJU0eSrX562i1LlwaghXZanyeHRSl+j0fQM0VMBmOW6m/255TTYFM6F+8hwjsLXw6WHB9d70Upfo9H0DAHREDSI+OqNVNfZ2JddSmhtCjX+g3p6ZL0arfQ1Gk3PETuLPoVb8KCaJT9tx1cq8YzQm7hdiT3pEiNEZKWIJInIHhH5tSkfJSLrRWSXiHwlIj7N2jwsIskisl9ELmwmn23KkkXkoa6ZkkajOWuInYXFVstk6x4OJBrZU8MG603crsSelX498Hul1FBgAnCPiAwFXgMeUkqNAD4H7gcwy64DhgGzgZdExGqmW3wRmAMMBa4362o0Gkel/0Rw8Wae526ibKkAeGl3zS6lXaWvlMpWSm0z35djJEUPAwYBq81q3wFXmu/nAR8opWqUUkeAZIy0iuOAZKXUYaVULfCBWVej0TgqTi4QM41zbdsYLBlUOPmDZ1BPj6pX0yGbvohEYeTL3Qjs4ZjSvhqIMN+HAenNmmWYsrbkGo3GkYm9EP/6fGZat1IbMLinR9PrsVvpi4gX8CnwG6VUGXA78EsR2Qp4A7WdMSARuVNEtojIlvz8/M7oUqPRnMnEzgTAXyrw6T+yhwfT+7FL6YuIM4bCf1cp9RmAUmqfUmqWUmos8D5wyKyeybFVP0C4KWtL3gKl1EKlVIJSKiE4OLij89FoNGcb3v0gZBQATv30Nl9XY4/3jgCvA3uVUv9sJu9j/rUAfwJeMYsWA9eJiKuIRAOxwCZgMxArItEi4oKx2bu4Myej0WjOUmJNJ78+Wul3NfbE058E3ATsEpEdpuyPGAr8HvP6M+BNAKXUHhH5CEjC8Py5RynVACAi9wLfAlbgDaXUns6aiEajOYtJuB3qj0LYmJ4eSa9HlFI9PYY2SUhIUFu2bOnpYWg0Gs1ZhYhsVUoltFamT+RqNBqNA6GVvkaj0TgQWulrNBqNA6GVvkaj0TgQWulrNBqNA6GVvkaj0TgQWulrNBqNA6GVvkaj0TgQZ/ThLBHJB1JPo4sgoKCThnO24IhzBsectyPOGRxz3h2dc6RSqtXgZWe00j9dRGRLW6fSeiuOOGdwzHk74pzBMefdmXPW5h2NRqNxILTS12g0Ggeityv9hT09gB7AEecMjjlvR5wzOOa8O23Ovdqmr9FoNJqW9PaVvkaj0WiaoZW+RqPROBC9UumLyGwR2S8iySLyUE+Pp6sQkQgRWSkiSSKyR0R+bcoDROQ7ETlo/vXv6bF2NiJiFZHtIrLEvI4WkY3md/6hmZKzVyEifiLyiYjsE5G9IjKxt3/XIvJb89/2bhF5X0TceuN3LSJviEieiOxuJmv1uxWD5835J4pIh9KN9TqlLyJW4EVgDjAUuF5EemvizXrg90qpocAE4B5zrg8BK5RSscAK87q38Wtgb7Prp4F/KaUGAsXAHT0yqq7lOeAbpVQcMApj/r32uxaRMOA+IEEpNRwjzep19M7v+i1g9nGytr7bORi5x2OBO4GXO3KjXqf0gXFAslLqsFKqFvgAmNfDY+oSlFLZSqlt5vtyDCUQhjHfRWa1RcBlPTLALkJEwoGLgdfMawGmA5+YVXrjnH2B84DXAZRStUqpEnr5d42Rx9tdRJwADyCbXvhdK6VWA0XHidv6bucBbyuDDYCfiITYe6/eqPTDgPRm1xmmrFcjIlHAaGAj0FcplW0W5QB9e2pcXcS/gQcAm3kdCJQoperN6974nUcD+cCbplnrNRHxpBd/10qpTOBZIA1D2ZcCW+n933UjbX23p6XjeqPSdzhExAv4FPiNUqqseZkyfHJ7jV+uiMwF8pRSW3t6LN2MEzAGeFkpNRqo5DhTTi/8rv0xVrXRQCjgyYkmEIegM7/b3qj0M4GIZtfhpqxXIiLOGAr/XaXUZ6Y4t/Hnnvk3r6fG1wVMAi4VkRQM0910DFu3n2kCgN75nWcAGUqpjeb1JxgPgd78XV8AHFFK5Sul6oDPML7/3v5dN9LWd3taOq43Kv3NQKy5w++CsfGzuIfH1CWYtuzXgb1KqX82K1oM3GK+vwX4srvH1lUopR5WSoUrpaIwvtsflFI3AiuBq8xqvWrOAEqpHCBdRAabohlAEr34u8Yw60wQEQ/z33rjnHv1d92Mtr7bxcDNphfPBKC0mRmofZRSve4FXAQcAA4Bj/T0eLpwnpMxfvIlAjvM10UYNu4VwEHgeyCgp8faRfOfBiwx3w8ANgHJwMeAa0+PrwvmGw9sMb/vLwD/3v5dA38B9gG7gXcA1974XQPvY+xb1GH8qrujre8WEAwPxUPALgzvJrvvpcMwaDQajQPRG807Go1Go2kDrfQ1Go3GgdBKX6PRaBwIrfQ1Go3GgdBKX6PRaBwIrfQ1Go3GgdBKX6PRaByI/wfXxFDdB34RhAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#!/usr/bin/env python\n",
    "# -*- coding:utf-8 -*-\n",
    "# Created on 2022-07-12 13:39\n",
    "# Author: FATE ZHOU\n",
    "\n",
    "from __future__ import division, print_function # Loading modules\n",
    "import time\n",
    "import datetime\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from PyEMD import EMD, EEMD, CEEMDAN # CEEMDAN # pip install EMD-signal\n",
    "from sampen import sampen2 # Sample Entropy\n",
    "from vmdpy import VMD # VMD\n",
    "# Sklearn\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, mean_absolute_percentage_error # R2 MSE MAE MAPE\n",
    "from sklearn.preprocessing import MinMaxScaler # Normalization\n",
    "# Keras\n",
    "from tensorflow.keras.models import Sequential\n",
    "from tensorflow.keras.layers import Dense, Activation, Dropout, LSTM, GRU\n",
    "from tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping\n",
    "\n",
    "# 1.Decomposition function\n",
    "# ========================================================================\n",
    "def ceemdan_decompose(series=None, trials=10, num_clusters = 3): # CEEMDAN Decompose\n",
    "    decom = CEEMDAN()\n",
    "    decom.trials = trials # Number of the white noise input\n",
    "    df_ceemdan = pd.DataFrame(decom(series.values).T)\n",
    "    df_ceemdan.columns = ['imf'+str(i) for i in range(len(df_ceemdan.columns))]\n",
    "    return df_ceemdan\n",
    "\n",
    "def sample_entropy(df_ceemdan=None, mm=1, r=0.1): # Sample Entropy Calculate; mm = 1 or 2; r = 0.1 or 0.2 \n",
    "    np_sampen = []\n",
    "    for i in range(len(df_ceemdan.columns)):\n",
    "        sample_entropy = sampen2(list(df_ceemdan['imf'+str(i)].values),mm=mm,r=r,normalize=True)\n",
    "        np_sampen.append(sample_entropy[1][1])\n",
    "    df_sampen = pd.DataFrame(np_sampen, index=['imf'+str(i) for i in range(len(df_ceemdan.columns))], columns=[CODE])\n",
    "    return df_sampen\n",
    "\n",
    "def kmeans_cluster(df_sampen=None, num_clusters=3): # K-Means Cluster by Sample Entropy\n",
    "    np_integrate_form = KMeans(n_clusters=num_clusters, random_state=9).fit_predict(df_sampen)\n",
    "    df_integrate_form = pd.DataFrame(np_integrate_form, index=['imf'+str(i) for i in range(len(df_sampen.index))], columns=['Cluster'])\n",
    "    return df_integrate_form\n",
    "\n",
    "def integrate_imfs(df_integrate_form=None, df_ceemdan=None): # Integrate IMFs and Residue to be 3 Co-IMFs\n",
    "    df_tmp = pd.DataFrame()\n",
    "    for i in range(df_integrate_form.values.max()+1):\n",
    "        df_tmp['imf'+str(i)] = df_ceemdan[df_integrate_form[(df_integrate_form['Cluster']==i)].index].sum(axis=1)\n",
    "    df_integrate_result = df_tmp.T # Use Sample Entropy sorting the Co-IMFs\n",
    "    df_integrate_result['sampen'] = sample_entropy(df_tmp).values\n",
    "    df_integrate_result.sort_values(by=['sampen'], ascending=False, inplace=True)\n",
    "    df_integrate_result.index = ['co-imf'+str(i) for i in range(df_integrate_form.values.max()+1)]\n",
    "    df_integrate_result = df_integrate_result.drop('sampen', axis=1, inplace=False)\n",
    "    return df_integrate_result.T\n",
    "\n",
    "def vmd_decompose(series=None, alpha=2000, tau=0, K=10, DC=0, init=1, tol=1e-7, draw=True): # VMD Decomposition\n",
    "    imfs_vmd, imfs_hat, omega = VMD(series, alpha, tau, K, DC, init, tol)  \n",
    "    df_vmd = pd.DataFrame(imfs_vmd.T)\n",
    "    df_vmd.columns = ['imf'+str(i) for i in range(K)]\n",
    "    return df_vmd \n",
    "\n",
    "# 2.Forecasting function\n",
    "# ========================================================================\n",
    "def GRU_model(trainset_shape):# Build GRU model\n",
    "    model = Sequential()\n",
    "    model.add(GRU(128, input_shape=(trainset_shape[1], trainset_shape[2]), activation='tanh', return_sequences=True))\n",
    "    model.add(Dropout(0.2))\n",
    "    model.add(GRU(64,activation='tanh',return_sequences=True))\n",
    "    model.add(Dropout(0.2))\n",
    "    model.add(GRU(32,activation='tanh',return_sequences=False))\n",
    "    model.add(Dropout(0.2))\n",
    "    model.add(Dense(1,activation='tanh'))\n",
    "    model.compile(loss='mse', optimizer='adam')\n",
    "    return model\n",
    "        \n",
    "def evaluation_model(y_test, y_pred): # Model evaluation function\n",
    "    y_test,y_pred = np.array(y_test).ravel(),np.array(y_pred).ravel()\n",
    "    r2 = r2_score(y_test, y_pred)\n",
    "    rmse = mean_squared_error(y_test, y_pred, squared=False) # MSE and MAE are different on different scales\n",
    "    mae = mean_absolute_error(y_test, y_pred)\n",
    "    mape = mean_absolute_percentage_error(y_test, y_pred)\n",
    "    df_evaluation = pd.DataFrame({'r2': r2, 'rmse': rmse, 'mae': mae, 'mape': mape}, index = range(1))\n",
    "    return df_evaluation\n",
    "\n",
    "def create_train_test_set(data=None, timestep=30, co_imf_predict_for_fitting=None): # Create training set and test set with normalization\n",
    "    if isinstance(data, pd.DataFrame): # Initialize DataFrame training set and test set\n",
    "        dataY = data['sum'].values.reshape(-1, 1)\n",
    "        dataX = data.drop('sum', axis=1, inplace=False)\n",
    "    else: # Initialize Series\n",
    "        dataY = data.values.reshape(-1, 1)\n",
    "        dataX = dataY\n",
    "\n",
    "    scalarX = MinMaxScaler(feature_range=(0,1)) # Normalize by sklearn\n",
    "    dataX = scalarX.fit_transform(dataX)\n",
    "    if co_imf_predict_for_fitting is not None: co_imf_predict_for_fitting = scalarX.transform(co_imf_predict_for_fitting)\n",
    "    \n",
    "    scalarY = MinMaxScaler(feature_range=(0,1))\n",
    "    dataY = scalarY.fit_transform(dataY)\n",
    "    \n",
    "    trainX, trainY = [], [] # Create training set and test set\n",
    "    for i in range(len(dataY)-timestep):\n",
    "        trainX.append(np.array(dataX[i:(i+timestep)]))\n",
    "        trainY.append(np.array(dataY[i+timestep]))\n",
    "        if co_imf_predict_for_fitting is not None: # When fitting, it uses today's forecasting result \n",
    "            if i<(len(dataY)-timestep-len(co_imf_predict_for_fitting)): trainX[i] = np.insert(trainX[i], timestep, dataX[i+timestep], 0)\n",
    "            else: trainX[i] = np.insert(trainX[i], timestep, co_imf_predict_for_fitting[i-(len(dataY)-timestep-len(co_imf_predict_for_fitting))], 0)\n",
    "    \n",
    "    return np.array(trainX), np.array(trainY), scalarY\n",
    "\n",
    "def GRU_predict(data=None, epochs=100, predict_duration=100, fitting=None): # GRU forecasting function\n",
    "    trainX,trainY,scalarY = create_train_test_set(data, co_imf_predict_for_fitting=fitting) # Get training and test X Y\n",
    "    x_train,x_test = trainX[:-predict_duration],trainX[-predict_duration:] # Split training and test set\n",
    "    y_train,y_test = trainY[:-predict_duration],trainY[-predict_duration:]\n",
    "    train_X = x_train.reshape((x_train.shape[0], x_train.shape[1], x_train.shape[2])) # Convert to tensor \n",
    "    test_X = x_test.reshape((x_test.shape[0], x_test.shape[1], x_test.shape[2])) # Convert to tensor \n",
    "\n",
    "    model = GRU_model(train_X.shape) # Build the model # Use model.summary() to show the model structure\n",
    "    patience = epochs//10 \n",
    "    EarlyStop = EarlyStopping(monitor='val_loss', patience=5*patience, verbose=0, mode='auto') # Early stop at small learning rate\n",
    "    Reduce = ReduceLROnPlateau(monitor='val_loss', patience=patience, verbose=0, mode='auto') # Adaptive learning rate\n",
    "    history = model.fit(train_X, y_train, epochs=epochs, batch_size=16, validation_split=0.1, verbose=0, shuffle=True, callbacks=[EarlyStop,Reduce]) # Train the model\n",
    "    \n",
    "    y_test_predict = model.predict(test_X) # Predict\n",
    "    df_gru_evaluation = evaluation_model(y_test, y_test_predict) # Evaluate model\n",
    "    y_test_predict = y_test_predict.ravel().reshape(-1,1) \n",
    "    y_test_predict_result = scalarY.inverse_transform(y_test_predict) # De-normalize \n",
    "    y_test_raw = scalarY.inverse_transform(y_test)    \n",
    "    df_predict_raw = pd.DataFrame({'raw': y_test_raw.ravel(), 'predict': y_test_predict_result.ravel()}, index=range(len(y_test_raw))) # Output\n",
    "    df_train_loss= pd.DataFrame({'loss': history.history['loss'], 'val_loss': history.history['val_loss']}, index=range(len(history.history['val_loss'])))\n",
    "    return df_predict_raw, df_gru_evaluation, df_train_loss\n",
    "\n",
    "# 3.Main function\n",
    "# ========================================================================\n",
    "if __name__ == '__main__':\n",
    "    start = time.time()\n",
    "    CODE, PATH = 'sh.000001', 'D:\\\\Stock-LSTM\\\\' # code such as 'sh.000001'\n",
    "\n",
    "    # 1.Load raw data\n",
    "    df_raw_data = pd.read_csv(PATH+CODE+'.csv', header=0, parse_dates=['date'], date_parser=lambda x: datetime.datetime.strptime(x, '%Y%m%d'))\n",
    "    series_close = pd.Series(df_raw_data['close'].values,index = df_raw_data['date'])\n",
    "\n",
    "    # 2.CEEMDAN decompose\n",
    "    df_ceemdan = ceemdan_decompose(series_close) \n",
    "    # df_ceemdan.plot(title='CEEMDAN Decomposition', subplots=True)\n",
    "\n",
    "    # 3.Sample Entropy Calculate\n",
    "    df_sampen = sample_entropy(df_ceemdan) \n",
    "    # df_sampen.plot(title='Sample Entropy')\n",
    "\n",
    "    # 4.K-Means Cluster by Sample Entropy\n",
    "    df_integrate_form = kmeans_cluster(df_sampen) \n",
    "    # print(df_integrate_form)\n",
    "\n",
    "    # 5.Integrate IMFs and Residue to be 3 Co-IMFs\n",
    "    df_integrate_result = integrate_imfs(df_integrate_form, df_ceemdan)\n",
    "    # df_integrate_result.plot(title='Integrated IMFs (Co-IMFs) of CEEMDAN', subplots=True)\n",
    "\n",
    "    # 6.Secondary Decompose the high-frequency Co-IMF0 by VMD\n",
    "    df_vmd_co_imf0 = vmd_decompose(df_integrate_result['co-imf0']) # vmd decomposition (The number of dataset must be even)\n",
    "    # df_vmd_co_imf0.plot(title='VMD Decomposition of Co-IMF0', subplots=True)\n",
    "\n",
    "    # 7.Predict Co-IMF0 by matrix-input GRU\n",
    "    time0 = time.time()\n",
    "    df_vmd_co_imf0['sum'] = df_integrate_result['co-imf0']\n",
    "    co_imf0_predict_raw, co_imf0_gru_evaluation, co_imf0_train_loss = GRU_predict(df_vmd_co_imf0)\n",
    "    print('======Co-IMF0 Predicting Finished======\\n', co_imf0_gru_evaluation)\n",
    "    time1 = time.time()\n",
    "    print('Running time: %.3fs'%(time1-time0))\n",
    "    # co_imf0_predict_raw.plot(title='Co-IMF0 Predicting Result')\n",
    "    # co_imf0_train_loss.plot(title='Co-IMF0 Training Loss')\n",
    "\n",
    "    # 8.Predict Co-IMF1 and Co-IMF2 by vector-input GRU\n",
    "    co_imf1_predict_raw, co_imf1_gru_evaluation, co_imf1_train_loss = GRU_predict(df_integrate_result['co-imf1'])\n",
    "    print('======Co-IMF1 Predicting Finished======\\n', co_imf1_gru_evaluation)\n",
    "    time2 = time.time()\n",
    "    print('Running time: %.3fs'%(time2-time1))\n",
    "    # co_imf1_predict_raw.plot(title='Co-IMF1 Predicting Result')\n",
    "    # co_imf1_train_loss.plot(title='Co-IMF1 Training Loss')\n",
    "\n",
    "    co_imf2_predict_raw, co_imf2_gru_evaluation, co_imf2_train_loss = GRU_predict(df_integrate_result['co-imf2'])\n",
    "    print('======Co-IMF2 Predicting Finished======\\n', co_imf2_gru_evaluation)\n",
    "    time3 = time.time()\n",
    "    print('Running time: %.3fs'%(time3-time2))\n",
    "    # co_imf2_predict_raw.plot(title='Co-IMF2 Predicting Result')\n",
    "    # co_imf2_train_loss.plot(title='Co-IMF2 Training Loss')\n",
    "\n",
    "    # 9. Add 3 result to get the final forecasting result (instead fitting method )\n",
    "    duration = 100\n",
    "    series_add_predict_result = co_imf0_predict_raw['predict']+co_imf1_predict_raw['predict']+co_imf2_predict_raw['predict']\n",
    "    df_add_predict_raw = pd.DataFrame({'predict': series_add_predict_result.values, 'raw': series_close[-duration:].values}, index=range(duration))\n",
    "    df_add_evaluation = evaluation_model(series_close[-duration:],series_add_predict_result)\n",
    "    print('======'+CODE+' Predicting Finished======\\n', df_add_evaluation)\n",
    "    end = time.time()\n",
    "    print('Total Running time: %.3fs'%(end-start))\n",
    "    df_add_predict_raw.plot(title=CODE+' Predicting Result')\n",
    "    # pd.DataFrame.to_csv(df_add_predict_raw, PATH+CODE+'_predict_output.csv')\n",
    "\n",
    "    # 10.Fit 3 result to get the final forecasting result (instead adding method )\n",
    "    \"\"\"\n",
    "    df_co_imf_predict_raw =  pd.DataFrame({'co-imf0': co_imf0_predict_raw['predict'], 'co-imf1': co_imf1_predict_raw['predict'], 'co-imf2': co_imf2_predict_raw['predict']}, index=range(len(co_imf0_predict_raw)))\n",
    "    df_fitting_set = df_integrate_result\n",
    "    df_fitting_set['sum'] = series_close.values\n",
    "    df_predict_raw, df_gru_evaluation, df_train_loss = GRU_predict(df_fitting_set, fitting=df_co_imf_predict_raw)\n",
    "    print('======'+CODE+' Predicting Finished======\\n', df_gru_evaluation)\n",
    "    end = time.time()\n",
    "    print('Running time: %.3fs'%(end-time3))\n",
    "    print('Total Running time: %.3fs'%(end-start))\n",
    "    df_predict_raw.plot(title=CODE+' Predicting Result')\n",
    "    df_train_loss.plot(title=CODE+' Training Loss')\n",
    "    # pd.DataFrame.to_csv(df_predict_raw, PATH+CODE+'_predict_output.csv')\n",
    "    \"\"\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}