Simple Autoencoders in Keras - 784-32-784-tf1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import os\n",
    "#os.environ['KERAS_BACKEND'] = 'plaidml.keras.backend' # https://github.com/plaidml/plaidml/blob/master/docs/index.md\n",
    "\n",
    "from tensorflow.keras.models import Model\n",
    "from tensorflow.keras.layers import Input, Dense\n",
    "from tensorflow.keras.datasets import mnist\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# keras version still end up import tf.keras version which is incompatible, import directly instead\n",
    "from tensorflow.keras.callbacks import TensorBoard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoding_dim = 32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_img = Input(shape=(784,))\n",
    "encoded = Dense(encoding_dim, activation='relu')(input_img) # apply to previous layer = input_img\n",
    "decoded = Dense(784, activation='sigmoid')(encoded)\n",
    "\n",
    "autoencoder = Model(input_img, decoded)\n",
    "autoencoder.run_eagerly = True # plaidml / keras not running eager mode by default\n",
    "\n",
    "encoder = Model(input_img, encoded)\n",
    "\n",
    "encoded_input = Input(shape=(encoding_dim,))\n",
    "\n",
    "decoder_layer = autoencoder.layers[-1] # extract last layer\n",
    "\n",
    "decoder = Model(encoded_input, decoder_layer(encoded_input))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.compile(optimizer='adam', loss='binary_crossentropy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "60000\n",
      "10000\n",
      "x_train shape = (60000, 784)\n",
      "x_test shape = (10000, 784)\n"
     ]
    }
   ],
   "source": [
    "# Load mnist dataset\n",
    "(x_train, _),(x_test, _) = mnist.load_data()\n",
    "\n",
    "print(len(x_train))\n",
    "print(len(x_test))\n",
    "\n",
    "x_train = x_train.astype('float32')/255\n",
    "x_test = x_test.astype('float32')/255\n",
    "x_train = x_train.reshape(len(x_train), (np.prod(x_train.shape[1:])))\n",
    "x_test = x_test.reshape(len(x_test), (np.prod(x_test.shape[1:])))\n",
    "\n",
    "print(f'x_train shape = {x_train.shape}')\n",
    "print(f'x_test shape = {x_test.shape}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1595258334\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 22:04 - loss: 0.6930WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (6.101998). Check your callbacks.\n",
      "WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.909997). Check your callbacks.\n",
      "  4/235 [..............................] - ETA: 10:59 - loss: 0.6902WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.457000). Check your callbacks.\n",
      "235/235 [==============================] - 52s 223ms/step - loss: 0.2780 - val_loss: 0.1894\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 31s 133ms/step - loss: 0.1721 - val_loss: 0.1556\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 33s 139ms/step - loss: 0.1271 - val_loss: 0.1197\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 31s 133ms/step - loss: 0.1168 - val_loss: 0.1116\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 22s 94ms/step - loss: 0.0928 - val_loss: 0.0917\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 21s 90ms/step - loss: 0.0927 - val_loss: 0.0916\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 20s 85ms/step - loss: 0.0927 - val_loss: 0.0916\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 21s 88ms/step - loss: 0.0927 - val_loss: 0.0916\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 22s 94ms/step - loss: 0.0927 - val_loss: 0.0916\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 20s 86ms/step - loss: 0.0927 - val_loss: 0.0916\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 19s 79ms/step - loss: 0.0927 - val_loss: 0.0915\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 19s 80ms/step - loss: 0.0927 - val_loss: 0.0916\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 19s 81ms/step - loss: 0.0926 - val_loss: 0.0915\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 18s 77ms/step - loss: 0.0926 - val_loss: 0.0915\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 18s 76ms/step - loss: 0.0926 - val_loss: 0.0915\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 16s 67ms/step - loss: 0.0926 - val_loss: 0.0915\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 16s 70ms/step - loss: 0.0926 - val_loss: 0.0915\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 15s 66ms/step - loss: 0.0926 - val_loss: 0.0916\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 20s 87ms/step - loss: 0.0926 - val_loss: 0.0915\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 16s 68ms/step - loss: 0.0926 - val_loss: 0.0915\n",
      "1595258334\n",
      "Wall time: 22min 12s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import os\n",
    "import time\n",
    "\n",
    "CURRENT_TIME = int(time.time())\n",
    "print(CURRENT_TIME)\n",
    "\n",
    "try:\n",
    "    os.makedirs(f'autoencoder_mnist/logs-{CURRENT_TIME}')\n",
    "except Exception:\n",
    "    pass\n",
    "\n",
    "# Run this in a terminal\n",
    "#tensorboard --logdir=autoencoder_mnist/logs\n",
    "\n",
    "CALLBACKS=[\n",
    "    TensorBoard(log_dir=f'autoencoder_mnist/logs-{CURRENT_TIME}')\n",
    "]\n",
    "\n",
    "# Fit data\n",
    "autoencoder.fit(x_train, x_train, epochs=50, batch_size=256, shuffle=True, validation_data=(x_test, x_test), callbacks=CALLBACKS)\n",
    "\n",
    "print(CURRENT_TIME)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoded_imgs = encoder.predict(x_test)\n",
    "decoded_imgs = decoder.predict(encoded_imgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbSklEQVR4nO3daZRUxRnG8RoBFURHGRmQsENQZJVdAwrKEWUTFJQwMQYEMWLEhcUoUQTUc0BRRATJCQqIiLIoCBKVAIpIPBCWsB8wgCO74OgoyMB0Pnh8favobnpm+nbfuf3/fXqKqumupOd291zrrUoLhUIGAAAAAAAA/nJOsicAAAAAAACAM3HTBgAAAAAAwIe4aQMAAAAAAOBD3LQBAAAAAADwIW7aAAAAAAAA+BA3bQAAAAAAAHyoZEEGp6WlcT54koRCobR4PA6vYVIdCYVC5ePxQLyOycO1GAhciwHAtRgIXIsBwLUYCFyLAcC1GAhhr0VW2gCJsyfZEwBgjOFaBPyCaxHwB65FwB/CXovctAEAAAAAAPAhbtoAAAAAAAD4EDdtAAAAAAAAfIibNgAAAAAAAD7ETRsAAAAAAAAf4qYNAAAAAACAD3HTBgAAAAAAwIdKJnsCSE2DBw+WXLp0aauvYcOGknv06BHxMSZNmiT5888/t/pmzJhR1CkCAAAAAJBUrLQBAAAAAADwIW7aAAAAAAAA+BA3bQAAAAAAAHyIPW2QMLNnz5Ycba8aLT8/P2LfgAEDJLdv397qW7FiheS9e/fGOkUkWZ06daz2tm3bJA8aNEjyhAkTEjanVHbBBRdIHjt2rGR97RljzNq1ayX37NnT6tuzZ49HswMAAEiOSy65RHLVqlVj+hn3O9FDDz0kedOmTZJ37NhhjduwYUNhpogAYaUNAAAAAACAD3HTBgAAAAAAwIcoj4JndDmUMbGXROmSmH/+85+Sa9asaY3r0qWL5Fq1all9WVlZkp999tmYnhfJd9VVV1ltXR6XnZ2d6OmkvMsuu0xy//79Jbtli02bNpXcuXNnq2/ixIkezQ5akyZNJM+bN8/qq169umfPe+ONN1rtrVu3Sv7qq688e16cnf6MNMaYBQsWSL7//vslT5482Rp3+vRpbycWQJmZmZLffvttyatWrbLGTZkyRfLu3bs9n9cv0tPTrfa1114recmSJZLz8vISNiegOOjUqZPkrl27Wn1t27aVXLt27Zgezy17qlatmuTzzjsv4s+VKFEipsdHcLHSBgAAAAAAwIe4aQMAAAAAAOBDlEchrpo1aya5e/fuEcdt3rxZsrvc8MiRI5Jzc3Mln3vuuda41atXS27UqJHVl5GREeOM4SeNGze22j/88IPk+fPnJ3o6Kad8+fJWe9q0aUmaCQqqQ4cOkqMtsY43twSnb9++knv16pWweeBn+rPvlVdeiTju5Zdfljx16lSr7/jx4/GfWMDoU2OMsb/T6FKkgwcPWuOSVRKlT/gzxn6v1+WtO3fu9H5ixcxFF11ktXXJff369SW7p5hSauZveluFgQMHStal4MYYU7p0aclpaWlFfl73lFQgVqy0AQAAAAAA8CFu2gAAAAAAAPgQN20AAAAAAAB8KKl72rhHQOs6wn379ll9J06ckDxz5kzJBw4csMZRj5tc+ohgt/ZT13zr/Rf2798f02M/8sgjVvvKK6+MOHbRokUxPSaST9eE62NojTFmxowZiZ5OynnggQckd+vWzepr0aJFgR9PHyVrjDHnnPPrfxvYsGGD5E8++aTAjw1byZK/foR37NgxKXNw98p4+OGHJV9wwQVWn96jCt7Q11/lypUjjps1a5Zk/f0KkV166aWSZ8+ebfWVK1dOst5L6C9/+Yv3E4tg+PDhkmvUqGH1DRgwQDLfm8+UlZUl+emnn7b6qlSpEvZn3L1vvvnmm/hPDHGj3x8HDRrk6XNt27ZNsv5bCPGjj1zX79XG2Hus6mPajTEmPz9f8uTJkyV/9tln1jg/vE+y0gYAAAAAAMCHuGkDAAAAAADgQ0ktjxozZozVrl69ekw/p5d1fv/991ZfIpedZWdnS3b/t6xZsyZh8/CThQsXStZL1YyxX6ujR48W+LHd42NLlSpV4MeA/1xxxRWS3XIKdwk64u+FF16QrJeJFtatt94asb1nzx7Jd9xxhzXOLbPB2bVr107y1VdfLdn9PPKSe/SxLlstU6aM1Ud5VPy5x7s//vjjMf2cLj0NhUJxnVNQNWnSRLK7xF4bOXJkAmZzpnr16lltXVI+f/58q4/P1jPpcpkXX3xRckZGhjUu0vUyYcIEq63LvQvznRexcUthdKmTLnFZsmSJNe6nn36SnJOTI9n9nNLfSz/88EOrb9OmTZL//e9/S163bp017vjx4xEfH7HT2ykYY19j+rum+zsRq5YtW0o+deqU1bd9+3bJK1eutPr079zJkycL9dyxYKUNAAAAAACAD3HTBgAAAAAAwIe4aQMAAAAAAOBDSd3TRh/xbYwxDRs2lLx161arr27dupKj1RW3atVK8ldffSU50hF94eg6tsOHD0vWx1m79u7da7VTdU8bTe9fUVhDhgyRXKdOnYjjdC1puDb8a+jQoZLd3xmuI28sXrxYsj6Su7D00aa5ublWX7Vq1STrY2e/+OILa1yJEiWKPI+gc+u59bHNu3btkvzMM88kbE633HJLwp4LZ2rQoIHVbtq0acSx+rvNBx984NmcgiIzM9Nq33bbbRHH3n333ZL190av6X1sPv7444jj3D1t3P0gYczgwYMl6yPcY+Xu03bTTTdJdo8N1/vfeLkHRlBF22emUaNGkvVRz67Vq1dL1n9X7t692xpXtWpVyXovU2Pisw8gzqTvBwwcOFCye41ddNFFYX/+66+/ttqffvqp5P/9739Wn/4bRO+t2KJFC2ucfk/o2LGj1bdhwwbJ+tjweGOlDQAAAAAAgA9x0wYAAAAAAMCHkloetXTp0qhtzT2q7RfucaONGzeWrJc5NW/ePOZ5nThxQvKOHTskuyVbeqmUXpqOouncubNkfXTmueeea407dOiQ5L/+9a9W348//ujR7FBU1atXt9rNmjWTrK83YzgaMV6uu+46q3355ZdL1st7Y13q6y7/1MuT9dGZxhhz/fXXS452HPGf//xnyZMmTYppHqlm+PDhVlsvEddL8d0StXjTn33u7xbLxRMrWsmOyy0jQHTPP/+81f7DH/4gWX+/NMaYd955JyFzcrVp00ZyhQoVrL7XX39d8htvvJGoKRUbunTXGGP69OkTdtzGjRut9sGDByW3b98+4uOnp6dL1qVXxhgzc+ZMyQcOHDj7ZFOc+/3/zTfflKzLoYyxy4OjlQxqbkmU5m5/gfh79dVXrbYua4t2fLe+b/Df//5X8mOPPWaN03/Xu6655hrJ+nvo1KlTrXH6/oJ+DzDGmIkTJ0qeO3eu5HiXyrLSBgAAAAAAwIe4aQMAAAAAAOBDSS2Piodjx45Z7WXLloUdF630Khq99NgtxdJLsWbPnl2ox8eZdLmMuyRS0/+fr1ixwtM5IX7ccgotkaduBJ0uQ3vrrbesvmjLTTV9mpde8vnUU09Z46KVI+rHuOeeeySXL1/eGjdmzBjJ559/vtX38ssvS87LyzvbtAOlR48ekt0TC3bu3Ck5kSet6TI3txxq+fLlkr/99ttETSllXXvttRH73FNpopUn4kyhUMhq69/1ffv2WX1engBUunRpq62X/t93332S3fn27dvXszkFgS53MMaYCy+8ULI+bcb9zqI/n37/+99LdksyatWqJblixYpW33vvvSf55ptvlnz06NGY5p4KypYtK9ndAkFvo3DkyBGr77nnnpPMVgn+4X6v06c29evXz+pLS0uTrP8ucEvnx44dK7mw2ylkZGRI1qeYjhgxwhqnt2lxSysThZU2AAAAAAAAPsRNGwAAAAAAAB/ipg0AAAAAAIAPFfs9bbyQmZkp+ZVXXpF8zjn2PS59HDV1qIX37rvvWu0bb7wx7Ljp06dbbff4WxQPDRo0iNin9zVB0ZQs+evbe6x72Lh7Q/Xq1UuyWzceK72nzbPPPit53Lhx1rgyZcpIdn8PFixYIHnXrl2Fmkdx1bNnT8n6/yNj7M8nr+k9krKysiSfPn3aGjd69GjJqbb/UKLoI0p1drk1/uvXr/dsTqmmU6dOVlsfp673cnL3YIiV3kelbdu2Vl+rVq3C/sycOXMK9Vyp6rzzzrPaek+gF154IeLP6eODX3vtNcn6vdoYY2rWrBnxMfReK17uh1ScdevWTfKjjz5q9eljuPWx98YYk5OT4+3EUCju+9iQIUMk6z1sjDHm66+/lqz3lv3iiy8K9dx6r5oqVapYffpvy8WLF0t297HV3PnOmDFDspd7+bHSBgAAAAAAwIe4aQMAAAAAAOBDlEeFMXDgQMn6WFr3ePHt27cnbE5Bc9lll0l2l3frJau6JEMvuzfGmNzcXI9mh3jTy7n79Olj9a1bt07yRx99lLA54Wf6qGj3iNjClkRFosucdImNMcY0b948rs9VXKWnp1vtSKUQxhS+9KIw9HHtutxu69at1rhly5YlbE6pKtZrJZG/H0E0fvx4q92uXTvJlSpVsvr00et66XzXrl0L9dz6MdyjvLUvv/xSsnvkNKLTx3W7dPmbW8IfSbNmzWJ+7tWrV0vmu2x40Uo/9ffG7OzsREwHRaRLlIw5s7RaO3XqlOSWLVtK7tGjhzXuiiuuCPvzx48ft9p169YNm42xv+dWqFAh4py0gwcPWu1ElYWz0gYAAAAAAMCHuGkDAAAAAADgQ5RHGWN+97vfWW13l/Jf6J3MjTFm06ZNns0p6ObOnSs5IyMj4rg33nhDcqqdGhMk7du3l1yuXDmrb8mSJZL1qQyIH/fkO00vPfWaXvLvzinaHEeMGCH5zjvvjPu8/MQ90eQ3v/mN5FmzZiV6OqJWrVph/53PwcSLVoYRj5OL8LO1a9da7YYNG0pu3Lix1XfTTTdJ1qeiHD582Bo3bdq0mJ5bn0ayYcOGiONWrVolme9IBeO+n+pSNl2C6JZg6BMwu3fvLtk9bUZfi25f//79JevXesuWLTHNPRW4pTCavt6efPJJq++9996TzIl5/vGvf/3LautSav03gjHGVK1aVfJLL70kOVqpqC63ckuxoolUEpWfn2+158+fL/mBBx6w+vbv3x/z8xUFK20AAAAAAAB8iJs2AAAAAAAAPsRNGwAAAAAAAB9iTxtjTMeOHa12qVKlJC9dulTy559/nrA5BZGuF27SpEnEccuXL5fs1qqieGrUqJFktyZ1zpw5iZ5OSrj33nslu7W5ydKlSxfJV111ldWn5+jOV+9pE3Tff/+91dY1+XpPDWPs/aGOHj0a13lkZmZa7Uj7C6xcuTKuz4vwWrduLbl3794Rx+Xk5EjmKNz4OnbsmGT3aHvdHjZsWJGfq2bNmpL1XmDG2O8JgwcPLvJzpaqPP/7YautrR+9b4+4zE2lfDffxBg4cKPn999+3+n77299K1vtj6M/tVFe+fHnJ7ncCvffbE088YfUNHz5c8uTJkyXrY9aNsfdN2blzp+TNmzdHnFO9evWstv67kPfb6NxjuPV+UBdffLHVp/eW1fvOfvPNN9a4vXv3Sta/E/pvDmOMadGiRYHnO2XKFKv92GOPSdb7VSUSK20AAAAAAAB8iJs2AAAAAAAAPpSy5VGlS5eWrI+OM8aYkydPStblOXl5ed5PLEDco7z10jJdgubSS39zc3PjPzEkRMWKFSW3adNG8vbt261x+hg9xI8uRUokvaTZGGOuvPJKyfo9IBr3mNxUeu91lxDrY3xvu+02q2/RokWSx40bV+Dnql+/vtXWJRnVq1e3+iKVBPil9C7o9OfpOedE/u9tH330USKmA4/pkg/32tPlV+57JWLnlpTefvvtknXZdnp6esTHmDBhgmS3LO7EiROS582bZ/Xp8o8OHTpIrlWrljUulY9xf+655yQ//PDDMf+cfn+87777wuZ40def3tqhV69ecX+uIHPLjfT1URjTp0+32tHKo3RJuv49e/31161x+kjxZGGlDQAAAAAAgA9x0wYAAAAAAMCHuGkDAAAAAADgQym7p82QIUMku0fPLlmyRPKqVasSNqegeeSRR6x28+bNw4579913rTbHfAfDn/70J8n6+OAPPvggCbNBojz++ONWWx97Gs3u3bsl33XXXVafPtYx1ej3Q/fo306dOkmeNWtWgR/7yJEjVlvvnXHppZfG9Bhu3Te8EenIdXcvgFdffTUR00Gc9ezZ02r/8Y9/lKz3XDDmzGNvER/6yG59vfXu3dsap685vfeQ3sPGNWrUKKtdt25dyV27dg37eMac+VmYSvS+JrNnz7b63nzzTcklS9p/ylapUkVytP2/4kHv4ad/Z/Sx48YYM3r0aE/nAWOGDh0quSB7Ct17772SC/M9KpFYaQMAAAAAAOBD3LQBAAAAAADwoZQpj9LLyI0x5m9/+5vk7777zuobOXJkQuYUdLEe0Xf//fdbbY75DoZq1aqF/fdjx44leCbw2uLFiyVffvnlhXqMLVu2SF65cmWR5xQU27Ztk6yPpDXGmMaNG0uuXbt2gR9bH2vrmjZtmtXOysoKO849ohzxUblyZavtlmj8Ijs722qvWbPGsznBOzfffHPEvvfff99q/+c///F6OilPl0rpXFju+6Qu99HlUe3atbPGlStXTrJ7RHnQ6SOW3fe1OnXqRPy5G264QXKpUqUkjxgxwhoXacuGwtLly02bNo3rYyO8fv36SdYlaW7JnLZ582arPW/evPhPzCOstAEAAAAAAPAhbtoAAAAAAAD4UKDLozIyMiS/9NJLVl+JEiUk66X9xhizevVqbycGi17+aYwxeXl5BX6MnJyciI+hl0emp6dHfIyLL77Yasda3qWXcA4bNszq+/HHH2N6jCDq3Llz2H9fuHBhgmeSmvRS3WgnKERblj9lyhTJlSpVijhOP35+fn6sU7R06dKlUD+XytavXx82x8OXX34Z07j69etb7U2bNsV1HqnqmmuusdqRrmH39EUUT+778A8//CD5+eefT/R04LG3335bsi6PuuOOO6xxevsAtm6IzdKlS8P+uy4nNsYujzp16pTk1157zRr397//XfKDDz5o9UUqW4U3WrRoYbX1e2PZsmUj/pzedkOfFmWMMT/99FOcZuc9VtoAAAAAAAD4EDdtAAAAAAAAfIibNgAAAAAAAD4UuD1t9F41S5YskVyjRg1r3K5duyTr47+ReBs3bizyY7zzzjtWe//+/ZIrVKgg2a0XjrcDBw5Y7aefftrT5/OT1q1bW+2KFSsmaSYwxphJkyZJHjNmTMRx+jjZaPvRxLpXTazjJk+eHNM4JIfeEylc+xfsYeMNvSef68iRI5LHjx+fiOnAA3pvBf09xRhjDh06JJkjvoNHf07qz+dbbrnFGvfkk09Kfuutt6y+HTt2eDS7YPrwww+ttv5+ro+I7t+/vzWudu3aktu2bRvTc2VnZxdihjgbd+/DCy+8MOw4vSeYMfa+UZ999ln8J5YgrLQBAAAAAADwIW7aAAAAAAAA+FDgyqNq1aoluWnTphHH6eOcdakU4sc9St1d9hlPPXv2LNTP6WP+opV1LFiwQPKaNWsijvv0008LNY8g6N69u9XWpYrr1q2T/MknnyRsTqls3rx5kocMGWL1lS9f3rPnPXz4sNXeunWr5HvuuUeyLmGE/4RCoahteKtDhw4R+/bu3Ss5JycnEdOBB3R5lHt9LVq0KOLP6ZKASy65RLL+vUDxsX79eslPPPGE1Td27FjJzzzzjNV35513Sj5+/LhHswsO/V3EGPvY9dtvvz3iz7Vr1y5i3+nTpyXra/bRRx8tzBQRhn6/Gzp0aEw/M3PmTKu9fPnyeE4paVhpAwAAAAAA4EPctAEAAAAAAPAhbtoAAAAAAAD4ULHf06ZatWpW2z3S7Rfung76mFt449Zbb7XauhaxVKlSMT1GvXr1JBfkuO6pU6dK3r17d8Rxc+fOlbxt27aYHx8/K1OmjOSOHTtGHDdnzhzJugYY3tmzZ4/kXr16WX3dunWTPGjQoLg+r3vM/cSJE+P6+EiM888/P2If+yd4Q38u6v35XCdOnJCcl5fn6ZyQHPpzMisry+p76KGHJG/evFnyXXfd5f3E4Knp06db7QEDBkh2v1OPHDlS8saNG72dWAC4n1sPPvig5LJly0pu1qyZNS4zM1Oy+/fEjBkzJI8YMSIOs4Qx9uuxZcsWydH+dtTXgH5tg4SVNgAAAAAAAD7ETRsAAAAAAAAfKvblUfoIWWOMqVq1athxK1assNocX5p4Y8aMKdLP9+7dO04zQbzopfnHjh2z+vQx6ePHj0/YnHAm95h13dYlpe77aZcuXSTr13PKlCnWuLS0NMl6KSuKrz59+ljtb7/9VvKoUaMSPZ2UkJ+fL3nNmjVWX/369SXv3LkzYXNCcvTr10/y3XffbfX94x//kMy1GCyHDx+22u3bt5fsluYMGzZMsltCh7M7ePCgZP1dRx+lbowxrVq1kvzUU09ZfYcOHfJodqnt+uuvl1y5cmXJ0f5212WjuoQ4SFhpAwAAAAAA4EPctAEAAAAAAPChtIKUCaWlpfmipqh169aSFy9ebPXpHae1Fi1aWG136bHfhUKhtLOPOju/vIYpam0oFGp29mFnx+uYPFyLgcC1eBYLFy602uPGjZO8bNmyRE8nrCBfi5UqVbLao0ePlrx27VrJATidLWWvRf1dVp8EZIxdwjpp0iSrT5cinzx50qPZFUyQr0W/cE/HvfrqqyW3bNlSchFKlFP2WgySIFyLGzZskNygQYOI48aOHStZlwsGQNhrkZU2AAAAAAAAPsRNGwAAAAAAAB/ipg0AAAAAAIAPFcsjv9u0aSM50h42xhiza9cuybm5uZ7OCQCAoNBHoCLx9u3bZ7X79u2bpJnAKytXrpSsj7gFwunRo4fV1vt+1K5dW3IR9rQBfKFcuXKS09J+3aLHPWL9xRdfTNic/ICVNgAAAAAAAD7ETRsAAAAAAAAfKpblUdHo5YI33HCD5KNHjyZjOgAAAABQaN99953VrlGjRpJmAnhr3LhxYfOoUaOscfv370/YnPyAlTYAAAAAAAA+xE0bAAAAAAAAH+KmDQAAAAAAgA+lhUKh2AenpcU+GHEVCoXSzj7q7HgNk2ptKBRqFo8H4nVMHq7FQOBaDACuxUDgWgwArsVA4FoMAK7FQAh7LbLSBgAAAAAAwIe4aQMAAAAAAOBDBT3y+4gxZo8XE0FU1eL4WLyGycPrWPzxGgYDr2Pxx2sYDLyOxR+vYTDwOhZ/vIbBEPZ1LNCeNgAAAAAAAEgMyqMAAAAAAAB8iJs2AAAAAAAAPsRNGwAAAAAAAB/ipg0AAAAAAIAPcdMGAAAAAADAh7hpAwAAAAAA4EPctAEAAAAAAPAhbtoAAAAAAAD4EDdtAAAAAAAAfOj/KDYux1kFDUwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x288 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABEUAAAHECAYAAADfz5RiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAasElEQVR4nO3ba4xdBdn34XsObaeblnactqRtbLEpAq0aEBUMFP3AwUgKSjCRIEI0IJoYSDERMHKISlXUGg1GIrF+ECk2BElEAiLUKLEIGCtYoFIbaOlBWuhhhpnpdNjvBzJ53yZP38eJT+/1wH1dX/dK/quzuvbs/etqR7vdDgAAAIBqOps+AQAAAIAmiCIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAldY/n4M7OznZnZ05H6enpSdmJiBgeHo6RkZGOtMEGTZw4sd1qtVK2urvH9dfrP7Zr166d7XZ7ZupoQ6ZOndqeMWNGytbAwEDKTkTEvn37YnBwsMS9OH369Pbs2bNTtl5//fWUnTEbNmwocy9OnDixPXny5JStvXv3puyMabfbJe7FadOmtWfNmpWyNTQ0lLIzZsuWLWXuxZ6envbUqVNTtvr7+1N2IiJGRkbiwIEDJe7F3t7e9ty5c1O2Nm/enLIzZu/evWXuxc7OznbWd4DM74tDQ0Oxf//+Evdiq9VqT58+PWVrZGQkZSfije8ZQ0ND/+U1HG8UiSlTpvzPnNV/Y/HixSk7ERHr1q1L22paq9WKJUuWpGzNnJn73r9y5coXUgcbNGPGjLjxxhtTth5//PGUnYiIX/7yl2lbTZs9e3b87Gc/S9nK/AAfEXHGGWeUuRcnT54cp556asrW/fffn7JTzaxZs2LFihUpWxs2bEjZGXP11VeXuRenTp0a5513XsrW2rVrU3YiIjZu3Ji21bS5c+emfQ5YtmxZys6YBx54oMy92N3dHUcddVTK1vHHH5+yE5F73zdt+vTpcfnll6dsvfTSSyk7ERG/+tWvDvma/z4DAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlNQ9noNHR0djz549h+tcDvL888+n7EREDA8Pp2017aijjoply5albB199NEpO2NWrlyZutekvr6+uOSSS1K2/vKXv6TsRER0dHSkbTVty5Yt8eUvfzll67bbbkvZqWhwcDDWrVuXspV5L1500UVpW01rtVrx/ve/P2Xr4YcfTtmpaP/+/bF169aUra9+9aspOxER1113XdpW01566aW4/vrrU7YWLVqUsjPmgQceSN1r0vTp02Pp0qUpW0ceeWTKTkTEU089lbbVtG3btsXXv/71lK2zzjorZSfijd8Th+JJEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoqXs8B/f29sbZZ599uM7lIKecckrKTkTEd7/73bStpm3dujVuvPHGlK358+en7FS0ffv2WL58ecpWR0dHyk41s2fPjmuvvTZla82aNSk7FS1YsCBuv/32lK3zzz8/ZSfijd8VVezevTvuueeelK1Jkyal7FTU0dERnZ05/9bX39+fshMRMTo6mrbVtKlTp8bpp5+esvXggw+m7FQ0ceLEmDdvXsrWypUrU3YiIl599dW0raaddNJJ8cQTT6Rsbdy4MWUnIuJjH/vYIV/zpAgAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACV1tNvtf/vgKVOmtN/1rncdxtP5vwYGBlJ2IiKef/75GBwc7EgbbNDkyZPbCxcuTNk6cOBAys6YZ5999sl2u/2+1NGGdHR0/Ps37n9o3bp1WVNx4YUXxt///vcS92Kr1Wofe+yxKVs9PT0pO2PWrl1b5l7s6elpz58/P2XrnHPOSdmJiLjjjjtix44dJe7FzPfTzZs3Z01FRMTb3/72Mvdi5nvqZz/72ZSdiIhvf/vb8eKLL5a5F7u6ulK2pk+fnrIzZteuXWXuxcz31N/+9rdZU/GFL3whNmzYUOZezNo6+eSTs6bi6aefjv7+/v/yGnpSBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABK6h7PwaOjo7Fnz57DdS4H+da3vpWyExGxbNmytK2m9fX1xac//emUrWeeeSZlZ8yzzz6butekefPmxXXXXZeytWnTppSdiIjh4eG0raYdddRRceWVV6ZsZd+La9euTd1r0syZM+Oyyy5L2XryySdTdiIiRkZG0raaNm/evLjmmmtStiq9x2UbHByMv/71r02fBv+BxYsXx+rVq1O2brrpppSdMXfddVfqXpPe/e53x3333ZeydfPNN6fsRETs3r07batpPT09sXDhwpStW2+9NWUnIuLiiy8+5GueFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkrrHc/CBAwfilVdeOVzncpA1a9ak7ERE7Nu3L22raVOmTIlTTz01ZWvbtm0pOxWNjIyk/XxffvnllJ2IiMHBwbStpg0ODsbTTz+dsjUwMJCyU9GuXbvijjvuSNvK0t/fn7bVtNHR0bQ/7969e1N2Kurp6YkFCxakbK1atSplJyLSPnf/b7B379546KGHUrbmzJmTslPRzp074/bbb0/Z+tvf/payExHx2muvpW01berUqbFkyZKUrUcffTRlJ+L//9nGkyIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUFL3eA6eMWNGXHLJJYfrXA4yMDCQshMR0W6307aatm3btvja176WsnXOOeek7FTU29sbF1xwQcrWN7/5zZSdiIihoaG0rab19PTEokWLUrbe8573pOyM+fGPf5y616Surq6YNm1aytZrr72WshMR0dlZ599MJk+eHMcff3zK1t13352yU1FHR0dMnDgxZWv27NkpOxER69evT9tq2uuvvx79/f0pWytWrEjZqairqyt6e3tTto499tiUnYiI5557Lm2raXv27In77rsvZesDH/hAyk7EG383D6XOpx4AAACA/4coAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJTUPZ6Dd+/eHffdd9/hOpeDvPbaayk7ERF79uxJ22paX19fXHTRRSlbq1evTtmpaPPmzXHVVVelbPX19aXsRES02+20raZNmDAhZs2albJ17rnnpuxU1NnZGa1WK2XrQx/6UMpORMQ999yTttW0bdu2xTe+8Y2UrZNOOillp6IpU6bEBz/4wZSt4447LmUnIuKxxx5L22ravn374g9/+EPKVtbOmCVLlqTuNWloaCieeeaZlK3t27en7EREjIyMpG01raurK972trelbK1ZsyZlJ+KN95hD8aQIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJTUPZ6Dh4eHY8OGDYfrXA6yaNGilJ2IiK6urrStpu3YsSN++MMfpmwdccQRKTsVTZ8+PZYuXZqyddppp6XsREQ89dRTaVtN+9e//hU/+tGPUra2bduWslPR/v37Y9OmTSlbf/rTn1J2IiL27duXttW0vr6+uPTSS1O2lixZkrIz5tZbb03da1JHR0f09PSkbGXei/39/WlbTZs2bVp85CMfSdm6/vrrU3Yq2rt3bzz44IMpW2effXbKTkTE448/nrbVtCOPPDLOPPPMlK3M31NDQ0OHfM2TIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlNQ9noMnTZoUCxcuPFzncpC+vr6UnYiI7u5x/Rje1Pr6+uLiiy9O2XrooYdSdiqaNm1aLF26NGVrdHQ0ZScioqurK22raT09PXHMMcekbE2ZMiVlZ8zq1atT95o0YcKEmD17dspW5nVcv3592lbTWq1WnHjiiSlbixYtStmpqLOzMyZNmpSytXbt2pSdiIiBgYG0rabt3bs37bPj5ZdfnrIz5pFHHknda9KcOXPihhtuSNlas2ZNyk5ERLvdTttq2iuvvBKrVq1K2brzzjtTdiIili1bdsjXPCkCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACV1j+fgxYsXxxNPPHG4zuUg9957b8pORMTVV1+dttW0zZs3xxe/+MWUrfe9730pOxXt2LEjbrnllpSt7du3p+xERGzZsiVtq2mdnZ3RarVStrLf41avXp2616R3vvOd8fDDD6dsZf5cr7nmmrStpo2MjMS2bdtSti688MKUnYrmzp0by5cvT9nq6elJ2YmIuO2229K2mjZ79uy49tprU7Yuv/zylJ2KWq1WvPe9703ZeuSRR1J2qlmwYEH85Cc/SdnKagsREcPDw4d8zZMiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFBS93gOfvbZZ+OUU045XOdykBdeeCFlJyJi586daVtN6+npiYULF6ZsHXPMMSk7Y5544onUvSb19vbGJz7xiZStI444ImUn4o33mCrmzp0by5cvT9l6xzvekbJT0ZYtW+JLX/pSytYPfvCDlJ2IiAMHDqRtNa2npyeOO+64lK2rrroqZWfMqlWrUveatHHjxjj//PNTthYtWpSyExHRbrfTtpq2cePGtM82Z511VsrOmKeffjp1r0mTJk1K+67R29ubshMR0dXVlbbVtJGRkdi+fXvK1osvvpiyExGxf//+Q77mSREAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKKl7PAdPmTIlTjvttMN1LgeZM2dOyk5ExCOPPJK21bRJkybFggULUrZ+8YtfpOyMufPOO1P3mvTSSy/FV77ylZStyZMnp+xEvPHnqmL9+vVx4oknpmzdddddKTtjTj755NS9Ju3evTvuvffelK3f/e53KTsREZdddlnaVtN27NgR3/nOd1K2urq6UnYqmjVrVlx55ZUpWzt37kzZiYhotVppW03r6OiICRMmpGx95jOfSdkZs3LlytS9Jm3fvj1uueWWlK3zzjsvZSci4v7770/balqr1YoTTjghZevxxx9P2Yl44z3mUDwpAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAldY/n4MHBwVi/fv3hOpeDzJs3L2UnIqK7e1w/hje1VqsVJ5xwQspWb29vyk5Fw8PDsWnTppSt6667LmUnIuL5559P22razJkz44orrkjZeu6551J2Kpo2bVp89KMfTdn68Ic/nLITEfH666+nbTVtaGgo/vGPf6RsXXDBBSk7FY2OjsaePXtStm644YaUnYiIrVu3pm01bf78+XHrrbembLVarZSdiubMmZN2j2R+X9y+fXvaVtP++c9/xqc+9amUrdNPPz1l57/jSREAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAErqHs/BXV1dMXXq1MN1LgdZvHhxyk5ExG9+85u0rabt3LkzfvrTn6ZsXXnllSk7Y2666abUvSYdffTR8b3vfS9la2hoKGUnImLChAlpW00bHR2NV199NWXr+9//fspORTNnzowrrrgiZevPf/5zyk5ExFNPPZW21bQZM2bEpZdemrK1b9++lJ2KXnnllfj5z3+esnXuueem7ERErFy5Mm2raUceeWSceeaZKVs333xzyk5FAwMD8dhjj6VsnXjiiSk7ERG///3v07aatmDBgli1alXK1qOPPpqyExFx9913H/I1T4oAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEnd4zl4aGgoNmzYcLjO5SCTJk1K2YmIGBgYSNtqWl9fX1x88cUpW8PDwyk7Fe3YsSNWrFiRsvXyyy+n7EREbN26NW2raf39/bF27dqUrRtuuCFlZ8yFF16YutekDRs2xBlnnJGy9fnPfz5lJyLixRdfTNtq2u7du+PXv/51yta6detSdirq7e2NCy64IGWr1Wql7ERErF69Om2raf39/fHHP/4xZWvevHkpOxW98MIL8bnPfS5l65Of/GTKTkTEk08+mbbVtI0bN8bHP/7xlK3Mz4z79+8/5GueFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkjra7fa/f3BHx8sR8cLhO53GzG+32zObPokMb+FrGOE6vhW4hm8NruObn2v41uA6vvm5hm8NruObn2v45nfIaziuKAIAAADwVuG/zwAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAl/R883kFsEm+ovgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd9RU1dXH8e2bYgMLRVABRQG7IsauSKxYwCxLoliiJsZlScxyrawYYxIxRWNd0VgT01BiX4pGLEDESERUEAGFKEVBQEUERY2ahPcPV7a/vfPMZRjmgcvM9/PXhnOeeea5d865d+46++w1li5dagAAAAAAACiX/1vVbwAAAAAAAAD/i4c2AAAAAAAAJcRDGwAAAAAAgBLioQ0AAAAAAEAJ8dAGAAAAAACghHhoAwAAAAAAUEKfX57Oa6yxBvXBV5GlS5euUY/X4RyuUguWLl3asR4vxHlcdRiLDYGx2AAYiw2BsdgAGIsNgbHYABiLDaHFschKG2DleXVVvwEAZsZYBMqCsQiUA2MRKIcWxyIPbQAAAAAAAEqIhzYAAAAAAAAlxEMbAAAAAACAEuKhDQAAAAAAQAnx0AYAAAAAAKCEeGgDAAAAAABQQjy0AQAAAAAAKKHPr+o3gMbyf//32XPAtdZaK7QdcsghHp9yyike9+nTJ/Rbe+21PZ43b57H48aNC/1eeeUVj+fPnx/aRo0a5fGcOXM8/s9//lP4/tUaa6xR08+hdnrMi9qWLl1asV9RG2qnx1/HuVk85vn4cz4AAFh+1d73oFyK7mU5j6gVK20AAAAAAABKiIc2AAAAAAAAJcRDGwAAAAAAgBJiTxuskJy3qfvRHHzwwaHtoosu8rhHjx4er7nmmqGf7pfRrl07j7fbbrvQ71//+pfHs2fPDm0LFy70eO7cuR7/+9///t8/ogLyTlvH5z73ufDvDTbYwOOdd945tHXv3t3jKVOmeDx16tTQ79133/U47z/EfkTVy+emY8eOHh955JEe77nnnqHfggULPB46dGhomzRpkseffPJJXd5ns8rzrc6VGhfNXToeWmNskMu/aunnQK/HZmbrrruux++//77HH374YejHnLn8ivb8UjoGOM7I8vz5+c9/9jXti1/8YmjT66nGzTbPrsxrTv5dek7atGnj8cYbbxz66R6fer9qFvfd1Lm42c4jlo2VNgAAAAAAACXEQxsAAAAAAIASIj0KK6Ro+d7mm28e/q3LPD/44AOPP/7449BPlwdqClRelviFL3zB45kzZ4a2l156qcXXwKqXPzO6bHT33XcPbV27dvV48eLFHk+bNq3ia7KkdPnouNJzYWbWv39/j8855xyPO3XqFPq98847Hmsam5nZ5MmT6/I+Eec8s5ha2LZtW49zyul7773nsaayffTRR6FfUbqGfk40jS5/ZlROR9W5vuh3MYarl6+LmgJ16qmnhrZdd93V40cffdTj+++/P/TTzwvnomXrrLNO+Lem8uo50WNpFsffP//5T4/zeKjluOe0LJ0H1ltvvdCm91maKrc8KeTNqlI6Tj3GSn4N/Xc+v3o90H75nrcRx7CeA/1uYRaPk36ei8aY/oxeS83MjjjiCI/POOOM0NazZ0+P9VqYU811jOX0qOHDh3t83XXXeZy3fchprGg+rLQBAAAAAAAoIR7aAAAAAAAAlFCp0qOKdklXusQ6L+VcVcsA83LJRlyO2JL8d2tKVLdu3UKbpixde+21Ho8YMSL002pP+vpbbbVV6KfpGl26dAltWnWq3ooquDTi0uKVuTN/Xpaq6RsvvPCCx5oqZdaYx31l0fOr6RNmZhdffLHHWg2haBl+v379QtuwYcM8fvvtt1fovTYjPdY5FWmLLbbwWNMzcr8JEyZ4rOeg1vGr1+rOnTuHNk3D0FQQM7M333zT46JqJ412/ay2spAu36/2GOT5uVevXh5/73vfC22aTqfvI6dHNdrxr1U+trvssovHF154YWjTlNGxY8d6fM8994R+Ov6KPhfVVpkqSm/VdOO+ffuGNr0f++tf/9ri+1vW725kelxzWuqGG27osX5X0TRhs5gSU+uY0uOfU260Opx+fnT7AbPmvj8qOo+aSqrVMQcPHhz66TUup2Lp6xedY02nbN++fWgbNGiQx/q96frrrw/9/v73v3usny0sm56nar/TlPE6yEobAAAAAACAEuKhDQAAAAAAQAnx0AYAAAAAAKCEVsqeNpprqTmFmoNvZrbnnnt6nPcvUa+++qrHudTzvHnzPNacv6Kyz7mMmvbV957LqGq/XNZRSzk2cj5pLnu50UYbefzKK6+EtiFDhnis+5NUW5L7rbfeCv/ea6+9PNb9NszMvvnNb3r83HPPeaz7KNRLUenBRi83Xm0+b6WfMTPbdNNNW4zN4udkzpw5Hjf6cV2ZdD8oHaNm8XwU5QFrbvj+++8f2s477zyPr7rqKo8XLlwY+pUxf7gM9Ljnvbr23Xdfj3fccUePX3zxxdBP9zjQPeFqPeb6njp06BDaevTo4fHUqVNDm5Y61etiI18jzarPp69lPs37XJx11lke5+uivr7ez1BK9jN6jHTPKLO4x1fe/0uPod7/vfzyyxX76X4lRZ+Laq2//vrh38ccc4zHOj+Yxf1vHn/8cY+beR7Wc6D7mPzgBz8I/bQEtN5TPvDAA6Gf7kkya9as0FbtXkF6PnTuNovXXd1bJ9/nNuL8WrTnU6Xvbbqnl5nZgAEDPD7//PM93mSTTUI/nWPz79LxrN8587ys7ynfv+o+UnqtznsTNdN9b6VrZtF+irqv2H777Rf6nXnmmR5vueWWoU3Hx6RJkzy+4447Qr9Ro0Z5/MYbb4Q2/Rzo69V7PmWlDQAAAAAAQAnx0AYAAAAAAKCEWiU9qqgkspb01fKJZmaHH364x3lZqqZV6dLTvKxXy41q6k4uJaxL13KZPi1Lqq+fl7tp+bU//vGPoU1TuGop41lmej7zclxdFvbYY4+FNk2XqqWM5NFHHx3+rcvwM13+VlRitRb5HDbTksV6yClkmuKRx/0jjzziMUv460fPweWXX+5xXhZcacl+UYnmnC6j6Rq6HPmKK64I/W6//XaP8zLwZqbnaocddghtWl5d++Uyw7Nnz/a4HkvldU7N57tSOVyzeF71fTR6WeGiMqIrWmJUy/6amX35y1/2ON+z6HEeMWKEx1zDPqOf7W233Ta0bbbZZhV/TtOgfvOb33is95Nmlc9x/v9q06W0X17237t3b4/zfdD48eM9XrRo0TLfXyPKx1jPr9575HtN/Tmd03Ts5X633HJLaJs2bZrHRXNyUXqUfnfRsf3RRx9VfL1GlD/bet3Rc5DvITWNWFPx27RpE/rpa9x2222h7corr/RYz0d+DU1H1LQ2s3ju9DUWLFgQ+jXyec1jUb/z6z1Fvgc67rjjPNZnCnmu1mcA+fOi469Pnz4eb7fddqGfplg9++yzoU3vZ6dPn+5xva+trLQBAAAAAAAoIR7aAAAAAAAAlBAPbQAAAAAAAEqoVfa0yTmxmi+mJcx03xezuOdJLlmnuZxaKjbnbGt+veah5hxCLdGdS/Fp2TbNn8u5+1oO9/nnnw9t+rc1Wo6w/j2aC20Wc0QXL14c2mrZt0CP+U9+8pPQpjmQes7MzH72s595zP4Y9bWiZb7bt28f2g4++GCPdU8qs1jmu9Z9LyqVDizaW6LR7bHHHh5r+dKi/Z80NzePKW3L+5jo3Lv11lt7fM0114R+WqLxggsuCG25vGIz0X3D+vfvH9q6devmsZaq1Px8s5gLr5/z5SkzrJ8NLT2u49cs5p/rvg1m8XPS6PvYqBXdtybT89azZ8/Q1rFjx4o/t2TJEo8ffPDBFX5PtVwLyk7/prxnn46BvD+G7hGke0jVco1cHlrGWPcPM4slcMeMGRPaRo8e7XG+324Wea/LoUOHeqzjKp8b/U4zc+bMFmOzuF+f7i9kFq9xuj9G0f42uU3392ym+dQsjsW8f4weC/1ukI/RU0895fG5557r8T777BP6TZw40eMpU6aEtkp7luTvP5XuQ7N6XyvKTL+/5+/Xuj/Uscce6/E222wT+uX7zf/KY0XLquszBLO4H5nev2y//fahn+5x1qtXr9D21ltveXzRRRe1+J7qgZU2AAAAAAAAJcRDGwAAAAAAgBJqlfSoTJd46XK+vMxMlxflpU2azqTLUvMyM10a3LVrV4/zsjgtyaVLhs1iWpWWbuzSpUvop8tS8xK5Rl6qqOczLxHWZfi1HgNdMnfzzTd7rCXzzOKyx9NPPz20Pf744y2+X6y4astAK13Kmkv2acpMXraoy8xrfX+VlvA30+dC01nMzO6++26P11lnnYo/p8vmJ0yY4PHkyZNDv4022sjjvMx144039lhT4/LvPfrooyu+vqZS1aNkdZnllF8tO6klLc3i0mBd6j1v3ryKr19rGoaWFNeUqMMOOyz002XheTw3axqGyvNOpZTEovlJPyMDBw4MbVoCPL+GlnnOpairVevnZ3W05pprhn/rcc/HYUVTVfLnoOj86xJ+vUfNaR2a+n/55ZeHNi0n3EzXQj1vu+66a2jL9yb/leetYcOGeXzppZd6vNNOO4V+l1xyicd5i4Yf//jHHg8aNMjjnFajmjmlO9P7jPzdTNPU9H4hHy89r6+//rrH9913X+in3zVqvf9o1ntPle9t9N7wwgsvDG2aEqXf/fL3dU1du+OOOzzWdHEzs/nz53usZdXz+zr11FM9zvdbWoY8/y26XUprXiNZaQMAAAAAAFBCPLQBAAAAAAAooZWSHqV0WVhe5qT/zkvQKi0ty8uQdBmbVnDKr1e0fHXGjBkeF1VT0ff70ksvVf36jaSoUlitNB2gT58+HmuKnJnZr3/9a4/vueeeur8PtKyWZZ66rLBv376hTZcNP/HEE6FNqwa1dhWORqPpLDfeeGNo03QmlZeBX3311R5fe+21Ffvp8tV33303tGmaqu6qf9BBB4V+mi51/PHHh7YbbrjB40Yf25reYhYrfeVqNroMXJfsa5qqWbwe1VrxR1NFjjrqKI9zOpxWaXjttddCW7MuCy9S7b2N0qo3hx56aGjTe5b8OfjVr37lcbWpatWmmzYK/fv0umUW59RceUjH6e9+9zuP83yor6/pjfl36TyX01t1PuzXr5/HuZqnVtzUyp5mzXOPmun46N69e2jTqoha6fbss88O/XSuVVoB0SymseVxtMkmm3jc6Ne0eshVR7/97W97nKvHXnfddR4XzVHaVqnKsVnzjpV6KKrGd8wxx3ic7/m0IphuPzJu3LjQ74wzzvBYU36Lzln+TOi9jX7/zOmxRddk/QySHgUAAAAAANBkeGgDAAAAAABQQjy0AQAAAAAAKKGVvqeNWp59ZirJuWn6GrWWWNOy4UVlvUeMGOGxlotb3t/XjDTnr1OnTqFt8ODBLf7Mk08+Gf6t5cCrzQnOuYacp+VXy94yum9NLrOpeeT3339/aMsl5av5XbkUX1HJx0amOfP9+/cPbXq8dM586KGHQj8di1rSNiuaa7WE6U033eTxIYccEvrpeevRo0do0z148j4pjUDPR+fOnUPblltu6XEuB3vvvfd6PHfuXI+XJ5+7WlqeU/O+814cY8eO9biofC0+Vel8FJUG79mzp8e9evWq+Nr5vmT06NHL/L1Znk+bqXRtvv7o36v70ZiZ7bvvvh6fd955Hj/zzDOhn95T6rUpl2DXOS+Xw9VzrmN95MiRoZ/+m305PqVzre6/ZRbvMX//+997/OCDD4Z+et50b6Pvfve7oZ/ugZSPv87X1e4v1Wx0zvvpT38a2gYMGOBx/m6Qv6tVo5afwbLp9aNbt26h7cgjj/RY9zQ0i+NU9xjSe0gzs4ULF7b4e/N9ib5e/h544IEHevyVr3ylxfee5e+cur9ga861rLQBAAAAAAAoIR7aAAAAAAAAlNAqTY+qxxKiotJatS7/1dJjmtYxb9680O/WW2/1uNo0jmamS0W7dOnicS5HrOkzumT4zjvvDP3y0lalyyr1/ObPHKUWl63WJfB6DjTFI5fZnDNnjsd5KXm1c4T+rvx+G30J/3/pMTCL85iW5DaLx2T+/PkeaxlNsziv1Xoc9Ry+9dZbHheNvdymKUONmB6l505LpJvFtKRZs2aFtlGjRnlc7+Xd+fN08MEHe6wpHjkFSsvhauojPlXrONLl3kcffbTHeo9iFsfOI488Etpy+elKqr2vasS5VeerSZMmhTadK3MJ4rXXXtvjk08+2ePTTjst9NPz89xzz3ms10Ezs7333tvjrbbaKrTpPY2e08suuyz0I+3mf+m8Nnv27NA2ZMgQj8eMGVPxNfR6+vOf/9xjTUnO8rl46aWXWnxPpPB/Ztttt/X4hBNOCG2aTrP99tuHtnztqqd8fiqlmuNT+r0v3/tvvvnmHufjqsdS57itt9469NNzrSmqvXv3Dv20JLd+/zSL6ab6ucpjT+fuCRMmhLbHHnusxX71xkobAAAAAACAEuKhDQAAAAAAQAmt9PSoeiz102VUeRlcLZVt8pLGAw44wGNd3j18+PDQ74UXXlju39tM8nI3XQo3dOhQj3fcccfQT5eRTpkyxWOt1pX7Zbp8WKs8FFXAQX3pssj999/fY02tMIvLkN94443QVu240n55iWqzjM08F/bp08fjvFzzo48+8lgrXrR2FTyda4uWgedUn0o7/zfKudX5Klc9UFpxxKz+abmVKr6ZmR177LEe6/vV66BZXDbcKOdnVcjjQ5dtH3TQQR7n9G5dBq4VFs1qW7ad59NGTwHQv2/q1Kmh7YorrvD48MMPD2269F/vOTRtyiym2Wu1xHxc99xzz4rvUceVplhNnz694s/gfxXNn3vttZfH7dq1C21f/epXPe7bt2/F19D7zVdeeaViP90SYPz48aFNx3MzzKc675111lke5+pCOu/l73CaTj1jxgyPa/l+uKw2vc/V+5ZGnyerpcch319quvx6660X2iqluH3rW98K/9axuf7663tcVPmpiH5GdOyZxbn2hz/8YWibNm2ax1SPAgAAAAAAaDI8tAEAAAAAACghHtoAAAAAAACU0Cot+b08Ku1jU1QmTOX8uLZt23qsezqYmW200UYev/jiix7fcMMNod+SJUs8boZc0+WV92a4/PLLPdZybDn3UMu7XXLJJR4XlfjONHdf85bJM209eSzqGNN9ovJYHDt2rMfVligt2g+lWcei5labxX1r8l5OWjpaywIXjY9K+8os6+d0f4dBgwZ5nOcHPW95r4G811Ejy59fnQ87dOgQ2nQ/MO1XtFeDfk7yOdB/H3XUUaFthx128FjH8MiRIyu+32bVGmV7N910U491f7hM9z2aOXPmCv/eZrtm6rnK8+ZDDz3k8ahRo0KbznN6/5HPvc7LumfimmuuGfp17NjRYy1Ja2b2wQcfeHzhhRd6TInv5bPhhhuGf5944oke6542eb8NvWfVz8jEiRNDv3Hjxnmsnw8zs/3228/jI444wuP77rsv9Lv00ks9XrRoUQt/RWPRPaAqXXPMivdf0zLsN910k8f5eqclxXVc5nlz8uTJHuf7rDZt2ng8Z84cjxcvXhz6Net9qc5JeY8w3bNol112CW1allvnRr0OmpntscceHuu5yN8ri/Yp0r2IdJ+dn/70p6HfXXfd5XEei3kfxtbCShsAAAAAAIAS4qENAAAAAABACa2W6VEql7DUJWj6M3lJm5Yg1tgsLln9wx/+4HFeMtdsy4aXl5bANDPr16+fx3o+dFmimdm5557rca3lY1szXaY1lr6vCkXLBVW1f19evrrZZpt53KNHD491abdZTI9iTNUuL6/v2rWrx/lcv/feex5Xu6xTz29eeqqvkdtOOOEEjwcMGOBxXi6un7OFCxeGNk2NXF3HWxH9m7QksFksk6npUGZm3/jGNzzWZcILFiwI/dZaay2Pdfn5/PnzQz89zieddFJo0xQBnbNzKdtayko3gmrLxFb7+c3jqH///h5rOkA+3prCk+faat4fPpOvR/q5z/ct9aZpNzk99C9/+YvHeo9UbXpr1ohzajU0hdvMbMstt/RY57s8FjX1ZfDgwR7ffffdoZ+m45x88smhbeDAgS2+D71emsX05dGjR4e2Rrxf0u0p9HtC0XUlf371O92hhx7qcU6j0nsafY08b2qaWy4DreXG77nnHo+vvvrq0C+nWjYLPa5632kWS2jrPJYVzV1a5vu8887z+Dvf+U7opyXj87jRey4dp2PGjAn9ypB+ykobAAAAAACAEuKhDQAAAAAAQAnx0AYAAAAAAKCESrunTVEOm+ajFeXiapuWTzQzO+ecczzWHEqzWMrxgQce8Li1c5gbgeb+nnnmmaFNy7HpORw+fHjop2XVqs3ZbY19ZjTfVfMh8/5IWuK2EXOMi+hxz3nf22+/vcd6/KZPnx76vfrqqx4XnTf9XXn/nGY77i3Jx6Rdu3Ye654mZmbdunXzuH379h7nUtGV9obKx1t/t+5lZGb2i1/8wmPNKc9jVvfFue2220JbziNvNJqvr+PBzOzhhx/2OJ8fLVnat29fj/P+RpqjP2PGDI/zWNS5TXP1zeI51nPXrLn6y6Pa65EeV71emsX9oPQ85bFx6623ely0D4T+rkbZp211piVuzeLeCrrXlJnZlVde6XGeEyrJ51j/Xe09dSPQv2/SpEmh7eabb/Z4t91283jKlCmhn+51qSWC8zHWeVj37zCL+7bp/m75+Ot3l3yPVct9T9nHuv5N+l0s70ej7zvvq6ZloXv37u1xPn5Kj0u+fmrpcd0TLv97iy228PjZZ58N/UaOHOlxs96vFn3Wat0LT+fGIUOGeKz7/ZnFc6rf2czMLrjgAo//9re/rfB7ak2stAEAAAAAACghHtoAAAAAAACUUGnTo7JaSjjrUrhcRk+Xu+nyRjOzX/7ylx7nZVQopkvqjz322NCm50NTIXSpqVnlJWk5/UOXlOaf0devVAY+v4YuqTQz22effTzu3Lmzx3Pnzg39tAxjLqGr76tsSyKLjkstS2a1vKVZLEGsSxMnT54c+lWb+lKUHtWaJd5XV5pCkVP6OnXq5PFpp53m8SWXXBL6aYlDPf55ifBWW23l8S233BLaNP20qPSxpgXl9KiyjZ16078vjwctO5mX6Wsq2jbbbONxXgau17iJEyd6nEura0pd0bVPx5/+TDMrmk+rpT+naYtmZj179mzxZ+bMmRP+/fLLL7f4nrDyVFviXe8/NI3ULJ7vP/3pT6FNS9RWK7+PSqkiZUwJqCe9N8zH8cYbb2wxzsek0jnNx/Tjjz/2+Mknnwxtun3AiSee6HFOzdH0N/28mFW+zy1S9jlB76E13eXOO+8M/RYsWOBxLtGtaUqa8qYp+2bxnlWPS0751fudnKalY32DDTbw+Oyzzw79NO2GlOL60TGh393zfYme38ceeyy03XvvvR6Xff5jpQ0AAAAAAEAJ8dAGAAAAAACghEqbHlXrEj5dqta9e3ePTz/99NBPl8VptSIzs6lTp67w+2gWOU3l61//use5KlelJcM77rhj6KfpAJo2oOk2ZnGp47hx40KbnkNNDTnssMNCvyOOOMJjTZkzixWPNI3g+eefD/0233xzj/MSzpkzZ3pc5hSPenzOdWmomdlBBx3ksS4bzrvq6xJf1C5Xt1uyZEnFvjr/6dLsp556KvT7xz/+4XGvXr08Puuss0K//fbbz+P11lsvtOkcoZ+zRYsWhX7HHHOMxzlltZnkeUKXx+dzrEvENe0pz8s6xnTJfh73+vo6d5nFKhz6+jmttNrUkEZWj79b093M4vyqn5FHH3009Gv0SmsrU1GaW26rNM/l8aw/p2mqu+66a+in18y8nL+We4n8mWzWlGL9W3MqRLWpEZVStfN50d/13nvvhTZNq581a5bHWrXKLM7XuQqktul7LzrXZad/U64KpYrO1YQJEzweOHCgx9ddd13o169fvxZ/r6aFm5l16NDB42pTX4tS+IsUbQOgbc1635zHgKYd7rzzzh7n86RVpgYPHhzaqq3AVwastAEAAAAAACghHtoAAAAAAACUEA9tAAAAAAAASqi0e9pUK+etaTm2q666ymMt2WwWy8pdf/31oa1ZcwVrkXMudd+LotxP3Wfm+9//fmjT/TK0/GEutac5onkPjMWLF3uspVPzviv6PvL71bxWLX+bS+FqydX8PnJubKPRY6a5w2axHLGOqfHjx4d+1ebnF+0TsDrlbLeWvN/JsGHDPM6lLjUvWM9T3pNJP79t27b1OJceLRrreq6mT5/ucf68TJs2zWPOZ8uKPvdF161qj6fuE/Dmm29W/N2a/5/Lhtda7hrF86nuQ6XHf+jQoaFftftyFH0mGH8tK9pvQulYycdSX2Prrbf2OO8BqPR+plZ5XDbyOS6ag+rxd1fat2t59hrScfrGG294/PLLL4d+up9K165dQ5vee+pnpFG+w9T6d+g5mTt3rscXXHBB6Hfttdd63K1bN4/XX3/90C+Xcq9ES3nfcMMNoa2W7wL5c1z2ctStRefa22+/PbT16dPHYz1e+XhfdtllHutejWar11zIShsAAAAAAIAS4qENAAAAAABACa2W6VG6BCqX/zrzzDM93nfffT3Oy+yuueYaj5u5vOyKysv1dAna3nvvHdp0aacuN8xpT23atKnqd+uStpz+1rFjxxZ/V1EZPl3aaBaX0A0ZMsTjhx9+OPSbN2+exzl1qtFp+tqgQYNCm6aeaaqFLulF/eSxqOUtdS40i6UuNe0il+uuJdUlp2lpudrTTjvNYy1XbbZ6LVEti3ofs6Jl4Do/arn2XLpdX6NZl3PXSu9ndt9994r99J4ll2avB8q2f2p5ymRXe5x0fOh8uM4664R+es/avXv30Pb00097XEuZarPGK/mt93Z675H/Nj1e1aYzVZv+Wyt97/k7zYYbbuhxnp+1jLimRzXC+awXPRaanm1mNmrUKI+/9rWvedyuXbvQTz8zOe1G7/lvuukmj7UUtVlt2wA0SprbijrssMM87t+/f2jTsanH+Lbbbgv9rr766hb7rW5YaQMAAAAAAFBCPLQBAAAAAAAoIYM3YqYAAAY2SURBVB7aAAAAAAAAlNBqv6fNtttuG9rOPvtsj3W/jWeffTb0u+WWWzwm/7N2+dhNmTLF4969e4e2ffbZx+Pjjz/e45122in069Kli8daWjjnIb799tstxmYx71Q/B2uvvXboN3v2bI9zKbkHH3zQY923oR6ldRuFlkbM50fzrUeMGOFxPcqXFpVRXVbfZqGf2ZNOOim0/fa3v/V4//339zjvq1DpuOZzPWPGDI91XzEzsyeeeMJjLVWM8tH5Nu9N9M4773g8a9Ysj99///3QT/dIyvn/zToWq6XlffU4msVjqXsz5PNULR3bzVQOul5qPUa6357eE+X9SnSO/dKXvhTa7r//fo8/+OCDiu+p2r2JGuF869+q+yLmv02vQfl6VGmfn3y9q8eeGLrvju5bk8tNt2/f3uO859z8+fM91tLWtZSXbgb5uIwePdpj3Tcqf6/88MMPW/wZM7Nhw4Z5PG3aNI/zHplYPjoOLr74Yo913JjFcTpmzBiPTz/99NCvUfYHYqUNAAAAAABACfHQBgAAAAAAoIRWy/QoTXfRNBuzuMxQl0PdfffdoR9L11qHLlXT9Bgzs+HDh7cY10qXw+alxdqmy+nyMnBdHptLZzbCkuHWUGnJdU4vGzlypMf33nuvx7WmyBSdj1pKUzcTLRFsFstb7rbbbh4feOCBod92223nsZZqv/7660M/TddolGWoZVTt57zauSu/nl5bX3/99dA2ceJEj7Vcuy4dX57fjU/ptWvjjTf2OKf86nVs0qRJHtc699U6n3J+V4yWE9ZSzznlpmge1c9M0Ws0WlnvamkaTNu2bUNbpRLBZvEeUON6HLs8pjTVqXPnzh7n9A+dB5YsWRLaNDW16D63mc59kXwcxo0b5/HkyZM97tSpU+inKaiakpbbOM61y+nAxx13nMeaupa/p+l9yimnnOJxPe5Dy3gdZKUNAAAAAABACfHQBgAAAAAAoIRWm/QoXQLap08fj3NVlHXXXddjrZ6iy+DMWMbWCPQcFi2FYyf9+tLj/uabb3qsFdnM4tJCPT95eWO93xOWTdNDtbqTxiiHeqdEFdFKNPmzoONWU+WmTp0a+rFcfPnocdW0px/96Eehn1YdeuGFFzzO1bvqoVnTalYGTXf585//7HHfvn1DPz3uTz/9dGirlNbdzOdKU510HiuqYJfvRfQ16lEhSuUUfk0HefXVVz3WVGOz4kpYeQsCFKu2kph+dzSrvgoblo9+r9cqaWZmBxxwQIv98ne9sWPHevzaa6+t8HvScZrPdb3nhFqw0gYAAAAAAKCEeGgDAAAAAABQQjy0AQAAAAAAKKHVZk8bzf884YQTPM55cJqP9u6773o8e/bs0I8cRWDF6diptZT3iv5eoJG15mc9v7aW737uuedC2/jx4z1uzb0fmpnuT/P444+Htkp7G9Xj88F8uvLMmzfP4/PPP99jLftsZrbDDjt4PGbMmNCm4xSfqrRXTdEeetXuF1arovLiuh9g0Ryq+6vonp1m8fuO7uPDeK4vjmfr0PGx6aabhjbd06Zt27Ye57L311xzjce1lvmutI9NGe9tWGkDAAAAAABQQjy0AQAAAAAAKKFSpUfpUilNhzIz69Chg8dbbLFFiz9jFlM0Bg4c6LEuSTUr57InAABWlaJyuFi1WKLfeHQ5/5w5c0Jb/jdaX73HmJYqrrV8sPbLaXHVpoEBZaSf2VmzZoU2TffT5wF33XVX6PfMM8/U9X2UHSttAAAAAAAASoiHNgAAAAAAACXEQxsAAAAAAIASKtWeNkXlg+fOnevxgAEDPM4l8LQcWK3lvwAAAAA0ptbeJ6oee2XonjaffPJJaGOfKzSKt99+O/y7a9euq+idlBsrbQAAAAAAAEqIhzYAAAAAAAAltLzpUQvM7NXWeCPLossAdYngokWLVsXbWdk2q+NrrbJzCM5jA+AcNgbO4+qPc9gYOI+rP87hSrAS0qE4j6s/zmFjaPE8rkFOJAAAAAAAQPmQHgUAAAAAAFBCPLQBAAAAAAAoIR7aAAAAAAAAlBAPbQAAAAAAAEqIhzYAAAAAAAAlxEMbAAAAAACAEuKhDQAAAAAAQAnx0AYAAAAAAKCEeGgDAAAAAABQQv8PlKnzmQkJCDwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x288 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 13.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n = 10\n",
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(x_test[i].reshape(28, 28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "n = 10\n",
    "plt.figure(figsize=(20, 8))\n",
    "for i in range(n):\n",
    "    # What does encoded images looks like?\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(encoded_imgs[i].reshape(16*1,2)) # encoded 64 # (8,8) # (4, 4*4)\n",
    "    #plt.imshow(decoded_imgs[i].reshape(4, 14*14).T) # input or decoded 784 # (28,28) # (4, 14*14)\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "n = 10\n",
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(decoded_imgs[i].reshape(28,28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0047, PSNR (dB) = 23.29\n",
      "MSE = 0.0129, PSNR (dB) = 18.91\n",
      "MSE = 0.0023, PSNR (dB) = 26.31\n",
      "MSE = 0.0090, PSNR (dB) = 20.45\n",
      "MSE = 0.0098, PSNR (dB) = 20.08\n",
      "MSE = 0.0012, PSNR (dB) = 29.07\n",
      "MSE = 0.0111, PSNR (dB) = 19.53\n",
      "MSE = 0.0193, PSNR (dB) = 17.15\n",
      "MSE = 0.0163, PSNR (dB) = 17.88\n",
      "MSE = 0.0081, PSNR (dB) = 20.90\n",
      "Wall time: 4 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://dsp.stackexchange.com/questions/38065/peak-signal-to-noise-ratio-psnr-in-python-for-an-image\n",
    "import math\n",
    "import numpy as np\n",
    "\n",
    "def mse(img1, img2):\n",
    "    return np.mean( (img1 - img2) ** 2 )\n",
    "\n",
    "def psnr(img1, img2):\n",
    "    mserr = mse(img1, img2)\n",
    "    PIXEL_MAX = 1.0\n",
    "    try:\n",
    "        return 20 * math.log10(PIXEL_MAX / math.sqrt(mserr))\n",
    "    except ZeroDivisionError:\n",
    "        return 'Same image'\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i]\n",
    "    contrast = decoded_imgs[i]\n",
    "    a = mse(original, contrast)\n",
    "    b = psnr(original, contrast) # can also use cv2.PSNR rather than defined function psnr\n",
    "    if b != 'Same image':\n",
    "        print(f'MSE = {a:.4f}, PSNR (dB) = {b:.2f}')\n",
    "    else:\n",
    "        print(f'MSE = {a:.4f}, PSNR too high')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9147\n",
      "0.8371\n",
      "0.9038\n",
      "0.9151\n",
      "0.8352\n",
      "0.9325\n",
      "0.8361\n",
      "0.6817\n",
      "0.8094\n",
      "0.9213\n",
      "Wall time: 201 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/imamun93/Image-Similarities-using-SSIM/blob/master/blog5_ssim.ipynb\n",
    "import numpy as np\n",
    "from tensorflow.image import ssim\n",
    "from tensorflow import convert_to_tensor\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = np.asarray(ssim(convert_to_tensor(original), convert_to_tensor(contrast), max_val=1)) # <class 'tensorflow.python.framework.ops.EagerTensor'>\n",
    "    print(f'{a[0]:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9862\n",
      "0.9703\n",
      "0.9950\n",
      "0.9726\n",
      "0.9636\n",
      "0.9981\n",
      "0.9744\n",
      "0.9266\n",
      "0.9579\n",
      "0.9778\n",
      "Wall time: 168 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/tensorflow/models/blob/master/research/compression/image_encoder/msssim.py\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from scipy.ndimage.filters import convolve\n",
    "\n",
    "def _FSpecialGauss(size, sigma):\n",
    "  \"\"\"Function to mimic the 'fspecial' gaussian MATLAB function.\"\"\"\n",
    "  radius = size // 2\n",
    "  offset = 0.0\n",
    "  start, stop = -radius, radius + 1\n",
    "  if size % 2 == 0:\n",
    "    offset = 0.5\n",
    "    stop -= 1\n",
    "  x, y = np.mgrid[offset + start:stop, offset + start:stop]\n",
    "  assert len(x) == size\n",
    "  g = np.exp(-((x**2 + y**2)/(2.0 * sigma**2)))\n",
    "  return g / g.sum()\n",
    "\n",
    "\n",
    "def _SSIMForMultiScale(img1, img2, max_val=255, filter_size=11,\n",
    "                       filter_sigma=1.5, k1=0.01, k2=0.03):\n",
    "  \"\"\"Return the Structural Similarity Map between `img1` and `img2`.\n",
    "  This function attempts to match the functionality of ssim_index_new.m by\n",
    "  Zhou Wang: http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "  Returns:\n",
    "    Pair containing the mean SSIM and contrast sensitivity between `img1` and\n",
    "    `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  img1 = img1.astype(np.float64)\n",
    "  img2 = img2.astype(np.float64)\n",
    "  _, height, width, _ = img1.shape\n",
    "\n",
    "  # Filter size can't be larger than height or width of images.\n",
    "  size = min(filter_size, height, width)\n",
    "\n",
    "  # Scale down sigma if a smaller filter size is used.\n",
    "  sigma = size * filter_sigma / filter_size if filter_size else 0\n",
    "\n",
    "  if filter_size:\n",
    "    window = np.reshape(_FSpecialGauss(size, sigma), (1, size, size, 1))\n",
    "    mu1 = signal.fftconvolve(img1, window, mode='valid')\n",
    "    mu2 = signal.fftconvolve(img2, window, mode='valid')\n",
    "    sigma11 = signal.fftconvolve(img1 * img1, window, mode='valid')\n",
    "    sigma22 = signal.fftconvolve(img2 * img2, window, mode='valid')\n",
    "    sigma12 = signal.fftconvolve(img1 * img2, window, mode='valid')\n",
    "  else:\n",
    "    # Empty blur kernel so no need to convolve.\n",
    "    mu1, mu2 = img1, img2\n",
    "    sigma11 = img1 * img1\n",
    "    sigma22 = img2 * img2\n",
    "    sigma12 = img1 * img2\n",
    "\n",
    "  mu11 = mu1 * mu1\n",
    "  mu22 = mu2 * mu2\n",
    "  mu12 = mu1 * mu2\n",
    "  sigma11 -= mu11\n",
    "  sigma22 -= mu22\n",
    "  sigma12 -= mu12\n",
    "\n",
    "  # Calculate intermediate values used by both ssim and cs_map.\n",
    "  c1 = (k1 * max_val) ** 2\n",
    "  c2 = (k2 * max_val) ** 2\n",
    "  v1 = 2.0 * sigma12 + c2\n",
    "  v2 = sigma11 + sigma22 + c2\n",
    "  ssim = np.mean((((2.0 * mu12 + c1) * v1) / ((mu11 + mu22 + c1) * v2)))\n",
    "  cs = np.mean(v1 / v2)\n",
    "  return ssim, cs\n",
    "\n",
    "\n",
    "def MultiScaleSSIM(img1, img2, max_val=255, filter_size=11, filter_sigma=1.5,\n",
    "                   k1=0.01, k2=0.03, weights=None):\n",
    "  \"\"\"Return the MS-SSIM score between `img1` and `img2`.\n",
    "  This function implements Multi-Scale Structural Similarity (MS-SSIM) Image\n",
    "  Quality Assessment according to Zhou Wang's paper, \"Multi-scale structural\n",
    "  similarity for image quality assessment\" (2003).\n",
    "  Link: https://ece.uwaterloo.ca/~z70wang/publications/msssim.pdf\n",
    "  Author's MATLAB implementation:\n",
    "  http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "    weights: List of weights for each level; if none, use five levels and the\n",
    "      weights from the original paper.\n",
    "  Returns:\n",
    "    MS-SSIM score between `img1` and `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  # Note: default weights don't sum to 1.0 but do match the paper / matlab code.\n",
    "  weights = np.array(weights if weights else\n",
    "                     [0.0448, 0.2856, 0.3001, 0.2363, 0.1333])\n",
    "  levels = weights.size\n",
    "  downsample_filter = np.ones((1, 2, 2, 1)) / 4.0\n",
    "  im1, im2 = [x.astype(np.float64) for x in [img1, img2]]\n",
    "  mssim = np.array([])\n",
    "  mcs = np.array([])\n",
    "  for _ in range(levels):\n",
    "    ssim, cs = _SSIMForMultiScale(\n",
    "        im1, im2, max_val=max_val, filter_size=filter_size,\n",
    "        filter_sigma=filter_sigma, k1=k1, k2=k2)\n",
    "    mssim = np.append(mssim, ssim)\n",
    "    mcs = np.append(mcs, cs)\n",
    "    filtered = [convolve(im, downsample_filter, mode='reflect')\n",
    "                for im in [im1, im2]]\n",
    "    im1, im2 = [x[:, ::2, ::2, :] for x in filtered]\n",
    "  return (np.prod(mcs[0:levels-1] ** weights[0:levels-1]) *\n",
    "          (mssim[levels-1] ** weights[levels-1]))\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = MultiScaleSSIM(original, contrast, max_val=1.0)\n",
    "    print(f'{a:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Simple Autoencoders in Keras - 784-32-784-tf1.h5\")"
   ]
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}