Deep Autoencoders in Keras - 784-256-128-64-32-64-128-256-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(256, activation='relu')(input_img) # apply to previous layer = input_img\n",
    "encoded = Dense(128, activation='relu')(encoded)\n",
    "encoded = Dense(64, activation='relu')(encoded)\n",
    "encoded = Dense(32, activation='relu')(encoded) # the bottleneck = encoding_dim\n",
    "\n",
    "decoded = Dense(64, activation='relu')(encoded)\n",
    "decoded = Dense(128, activation='relu')(decoded)\n",
    "decoded = Dense(256, activation='relu')(decoded)\n",
    "decoded = Dense(784, activation='sigmoid')(decoded)\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[-4](encoded_input)\n",
    "decoder_layer = autoencoder.layers[-3](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-2](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-1](decoder_layer) # extract last layer\n",
    "\n",
    "decoder = Model(encoded_input, decoder_layer)"
   ]
  },
  {
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1595258015\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 19:06 - loss: 0.6924WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (5.220998). Check your callbacks.\n",
      "  3/235 [..............................] - ETA: 12:46 - loss: 0.6908WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.625999). Check your callbacks.\n",
      "235/235 [==============================] - 90s 382ms/step - loss: 0.2408 - val_loss: 0.1632\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 73s 313ms/step - loss: 0.1485 - val_loss: 0.1349\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 74s 314ms/step - loss: 0.1283 - val_loss: 0.1206\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 70s 298ms/step - loss: 0.1184 - val_loss: 0.1138\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 81s 345ms/step - loss: 0.1131 - val_loss: 0.1105\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 92s 390ms/step - loss: 0.1092 - val_loss: 0.1066\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 87s 371ms/step - loss: 0.1057 - val_loss: 0.1025\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 84s 359ms/step - loss: 0.1026 - val_loss: 0.1008\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 87s 372ms/step - loss: 0.1002 - val_loss: 0.0989\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 85s 363ms/step - loss: 0.0984 - val_loss: 0.0971\n",
      "Epoch 11/50\n",
      "235/235 [==============================] - 91s 386ms/step - loss: 0.0968 - val_loss: 0.0952\n",
      "Epoch 12/50\n",
      "235/235 [==============================] - 94s 400ms/step - loss: 0.0954 - val_loss: 0.0939\n",
      "Epoch 13/50\n",
      "235/235 [==============================] - 79s 335ms/step - loss: 0.0942 - val_loss: 0.0929\n",
      "Epoch 14/50\n",
      "235/235 [==============================] - 74s 314ms/step - loss: 0.0931 - val_loss: 0.0930\n",
      "Epoch 15/50\n",
      "235/235 [==============================] - 68s 290ms/step - loss: 0.0921 - val_loss: 0.0914\n",
      "Epoch 16/50\n",
      "235/235 [==============================] - 66s 279ms/step - loss: 0.0913 - val_loss: 0.0905\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 63s 269ms/step - loss: 0.0905 - val_loss: 0.0899\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 65s 278ms/step - loss: 0.0898 - val_loss: 0.0890\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 56s 239ms/step - loss: 0.0892 - val_loss: 0.0891\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 57s 244ms/step - loss: 0.0886 - val_loss: 0.0883\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 55s 234ms/step - loss: 0.0882 - val_loss: 0.0875\n",
      "Epoch 22/50\n",
      "235/235 [==============================] - 65s 276ms/step - loss: 0.0876 - val_loss: 0.0875\n",
      "Epoch 23/50\n",
      "235/235 [==============================] - 57s 244ms/step - loss: 0.0873 - val_loss: 0.0871\n",
      "Epoch 24/50\n",
      "235/235 [==============================] - 60s 254ms/step - loss: 0.0869 - val_loss: 0.0865\n",
      "Epoch 25/50\n",
      "235/235 [==============================] - 71s 304ms/step - loss: 0.0866 - val_loss: 0.0869\n",
      "Epoch 26/50\n",
      "235/235 [==============================] - 57s 242ms/step - loss: 0.0862 - val_loss: 0.0859\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 58s 249ms/step - loss: 0.0859 - val_loss: 0.0855\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 61s 260ms/step - loss: 0.0855 - val_loss: 0.0855\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 62s 265ms/step - loss: 0.0852 - val_loss: 0.0851\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 59s 250ms/step - loss: 0.0848 - val_loss: 0.0849\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 50s 212ms/step - loss: 0.0846 - val_loss: 0.0847\n",
      "Epoch 32/50\n",
      "235/235 [==============================] - 58s 248ms/step - loss: 0.0842 - val_loss: 0.0842\n",
      "Epoch 33/50\n",
      "235/235 [==============================] - 60s 256ms/step - loss: 0.0839 - val_loss: 0.0839\n",
      "Epoch 34/50\n",
      "235/235 [==============================] - 52s 223ms/step - loss: 0.0837 - val_loss: 0.0837\n",
      "Epoch 35/50\n",
      "235/235 [==============================] - 56s 240ms/step - loss: 0.0835 - val_loss: 0.0838\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 63s 267ms/step - loss: 0.0832 - val_loss: 0.0840\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 49s 210ms/step - loss: 0.0830 - val_loss: 0.0833\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 59s 252ms/step - loss: 0.0828 - val_loss: 0.0833\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 49s 208ms/step - loss: 0.0826 - val_loss: 0.0830\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 58s 248ms/step - loss: 0.0824 - val_loss: 0.0833\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 59s 250ms/step - loss: 0.0823 - val_loss: 0.0826\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 50s 214ms/step - loss: 0.0821 - val_loss: 0.0822\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 39s 167ms/step - loss: 0.0820 - val_loss: 0.0826\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 39s 167ms/step - loss: 0.0819 - val_loss: 0.0826\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 39s 167ms/step - loss: 0.0817 - val_loss: 0.0821\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 42s 177ms/step - loss: 0.0816 - val_loss: 0.0821\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 45s 191ms/step - loss: 0.0816 - val_loss: 0.0819\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 37s 156ms/step - loss: 0.0814 - val_loss: 0.0819\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 34s 143ms/step - loss: 0.0813 - val_loss: 0.0817\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 10s 43ms/step - loss: 0.0812 - val_loss: 0.0819\n",
      "1595258015\n",
      "Wall time: 52min 24s\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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoded_imgs = encoder.predict(x_test)\n",
    "decoded_imgs = decoder.predict(encoded_imgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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+WH4yJAAAZUklEQVR4nO3bf6zWdf3/8efhHDgcDtcBDqAgSaXoXL+kcIEbDrQJLdPp6A+zXDNWSWi69UeR0XStqOZwuZVrYFiydGxmrrmAxRL7MbeoZKPmj035DAhEjxA/DnD4cX3+aGf7ssU3zyfP8wqet9u/5+0eb3zvda7Lu2/ams1mAAAAAFQzotU3AAAAANAKoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJHUO5uLu7uzl+/PjhupfTnDx5MmUnIuLAgQPR39/fljbYQuPGjWuef/75KVsHDx5M2Rm0Z8+eN5rN5uTU0Rbp7OxsjhkzJmXr2LFjKTsREQMDA3HixIkSZ7Grq6vZaDRStiZPzj0Wf/vb38qcxe7u7mZvb2/K1sDAQMpOxD8/F48cOVLiLI4aNarZ1dWVsnXq1KmUnUGHDh0qcxZHjx7dHDt2bMpW5nM8fPhwHDt2rMRZ7OzsbHZ3d6dsZZ/Ff/zjH2XO4oQJE5oXXHBBytYrr7ySshMRcfz48TLfUUeNGtUcPXp0ylZWW4iIePPNN+PQoUP/8hkOKYqMHz8+br/99rfnrv6NQ4cOpexERPzkJz9J22q1888/Px588MGUrWeffTZlZ9CKFSv+J3WwhcaMGRPz589P2dq+fXvKTkTEiy++mLbVao1GIz7xiU+kbH3hC19I2Rk0c+bMMmext7c37rrrrpStnTt3puxERDz22GNpW63W1dUVc+bMSdnKjMwREZs3by5zFseOHRvXXXddylZ/f3/KTkTEpk2b0rZarbu7OxYsWJCydfTo0ZSdQU899VSZs3jBBRfEz372s5Stm2++OWUnIvf7cKuNHj06Zs+enbJ1/fXXp+xERNx///1n/Jm/PgMAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACU1DGUiw8ePBjPPvvscN3LaX7961+n7FQzcuTImDZtWsrWihUrUnYqOnHiROzfvz9l6/nnn0/ZqaazszMuvvjilK2ZM2em7FR06tSp6O/vT9n6/ve/n7JTTbPZjJMnT6Zsbd68OWWnookTJ8Ztt92WsnX11Ven7FTTaDRi3rx5KVtf/OIXU3YqOnDgQPzmN79J2XrhhRdSdqoZMWJEjB49OmXrrrvuStn5d7wpAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAltTWbzbd+cVvbW7/4LNNsNttafQ8ZzuVnGBF/ajabV7T6JjKcy8/RWTwnOIvnAGfxnOAsngOcxXOCs3gOcBbPfmd6ht4UAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgpI6hXHzeeefFpz71qeG6l9M0Go2UnYiIVatWpW212mWXXRaPPPJIytaPfvSjlJ1Ba9asSd1rpZEjR8akSZNStubPn5+yExGxYcOGtK1WmzVrVmzZsiVla9u2bSk7g97//ven7rVSb29vLFy4MGXroosuStmJiHj44YfTtlqts7MzLrzwwpStefPmpewMqvQce3p6Ys6cOSlb06dPT9mJiPjFL36RttVqjUYjZs+enbJ1/PjxlJ1BmzdvTt1rpSlTpsTixYtTtu68886UnYiIBQsWpG212oQJE9L+vB/+8IdTdiIiHnjggTP+zJsiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFCSKAIAAACUJIoAAAAAJYkiAAAAQEmiCAAAAFBSx1Aubm9vj+7u7uG6l9Ncd911KTsREevWrUvbarVdu3bF17/+9ZStSy65JGWnog984AOxZcuWlK3169en7ERE2p/pv8G+ffvSfve88sorKTsVjR07NubNm5eylfX5GxHR1dWVttVq73vf+9J+97z88sspO4Mefvjh1L1WmjRpUixevDhl69SpUyk7ERGbNm1K22q1tra26OgY0n+a/J/19vam7FTUaDTiqquuStn64Q9/mLITEfH666+nbbVaV1dXvPe9703ZGjlyZMpOxD9/x5yJN0UAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoKSOoVzc1tYWnZ2dw3Uvpzl16lTKTjWNRiPmzZuXsvW73/0uZaeiN954I1atWpWytWvXrpSdiIhDhw6lbbVaX19frF27NmXrk5/8ZMpORf39/bFly5aUrc997nMpOxGR9ln/3+DNN9+Mxx9/PGXryJEjKTsVHT9+PPbs2ZOydeedd6bsRETcf//9aVutNn78+LjhhhtSthqNRsrOoHXr1qXutVJ/f39s3bo1Zeuyyy5L2YmIGD16dNpWqx09ejReeumllK2FCxem7EREdHScOX14UwQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAASuoYysVjxoyJWbNmDde9nObKK69M2anmyJEj8de//jVla8OGDSk7Fe3duzd+8IMfpGxt3bo1Zaeanp6eWLBgQcrWLbfckrJTUU9PTyxcuDBla/bs2Sk71fT19cUjjzySsuVzcfiMGzcuPvaxj6VsjRjh/ykOh4MHD8YzzzyTsrVu3bqUnYrOO++8uOOOO1K2uru7U3aq6ezsjIsvvjhl69Zbb03Z+Xf8VgcAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAktqazeZbv7it7a1f/B+69dZbs6bi6aefjr6+vra0wRbKfIa//e1vs6YiIuKqq676U7PZvCJ1tEWcxbNf5jN8/PHHs6YiIuLmm292FofB5z//+ayp+PnPfx6vv/66s/g2W7RoUdZUREQ88cQTzuIwWLZsWdZUrFmzJnbv3u0svs06OzuzpiIi4tixY87iMHjhhReypmLRokWxbds2Z/FttnTp0qypWLduXezdu/dfPkNvigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASR1DuXjcuHExd+7c4bqX03z6059O2YmIeO6559K2Wq3RaMScOXNStk6cOJGyU9GkSZPipptuStn6yEc+krITEfH73/8+bavVGo1GzJ49O2Xrj3/8Y8pORaNHj44ZM2akbPX29qbsRER0dAzp68FZbcqUKXHbbbelbF1zzTUpO4OeeOKJ1L1WyvyO+u1vfztlJyJi48aNaVut1tnZGRdeeGHK1vz581N2Bq1evTp1r5Xa29tj3LhxKVuvvvpqyk5ExMDAQNpWq02cODFuuOGGlK13vOMdKTsREaNGjTrjz7wpAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASR1DubizszMuvfTS4bqX06xatSplJyLijTfeSNtqtalTp8ayZctStubPn5+yU1F3d3dceeWVKVuTJ09O2YmI6OgY0q+ks1qz2YyjR4+mbI0dOzZlp6Lx48fHxz/+8ZSt7du3p+xERBw7dixtq9XGjBkTl19+ecrWJZdckrJT0YgRI6K7uztla/369Sk7EREHDhxI22q1KVOmxFe/+tWUrblz56bsDFq9enXqXiu9+93vjgcffDBlq6urK2UnIqKtrS1tq9UGBgZi586dKVv/LZ+L3hQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJI6hnLxiBEjorOzc7ju5TSNRiNlJ+Kff64q9u/fH08++WTK1ne+852UnYpee+21+N73vpeyNX/+/JSdiIh9+/albbXatGnTYsWKFSlbS5YsSdmpaOrUqbF8+fKUraVLl6bsRNT6XPz73/8e9957b8rWZz7zmZSdinp7e+OWW25J2dq6dWvKTkREf39/2larHTlyJLZt25ayNXHixJSdirq6uuI973lPytZjjz2WshMRcfjw4bStVuvu7o5Zs2albO3YsSNlJyJiYGDgjD+r860HAAAA4P8higAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAldQzl4mnTpsWKFSuG615Oc/fdd6fsRES0t7enbbXayZMn49ChQylb9913X8rOoI0bN6butVKz2YwTJ06kbO3evTtlJyJiYGAgbavVjh8/Hq+99lrK1k9/+tOUnUEf+tCHUvdaqa+vLx599NGUrQkTJqTsRNT6XOzp6Ylrr702ZeuZZ55J2anowIEDsWHDhpStw4cPp+xERBw9ejRtq9Xa29ujp6cnZevGG29M2ako8/vNBz/4wZSdiIgxY8akbbVad3d3zJ49O2Vrz549KTsREU8//fQZf+ZNEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoqa3ZbL71i9va3vrF/6GVK1dmTcXKlStjx44dbWmDLZT5DL/1rW9lTUVExD333POnZrN5Repoi2Q+x1/96ldZU/GlL30pXnrpJWfxbfbd7343ayoiIr7yla84i8Pgy1/+ctZUrF27Nvbs2eMsvs2OHDmSNRUREV1dXc7iMFizZk3WVNx3333x6quvOotvs5kzZ2ZNRUTE888/7ywOg3379mVNxdVXXx1/+ctfnMW32VBaxH/qiiuuiC1btvzLZ+hNEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoSRQBAAAAShJFAAAAgJJEEQAAAKAkUQQAAAAoqWMoF0+aNCkWLVo0XPdymosuuihlJyKis7MzbavV3vnOd8by5ctTtlavXp2yU9G0adPijjvuSNlau3Ztyk5ERF9fX9pWq11++eWxcePGlK0f//jHKTsVTZ8+PZYtW5aytWnTppSdiIiBgYG0rVabPn163HPPPSlbWb+3GV4vvvhi2tbRo0fTtlrtXe96V9x7770pW6dOnUrZGfTZz342da+VLr300njooYdStmbMmJGyExGxf//+tK1WmzZtWixdujRla/v27Sk7Ef//7zbeFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoCRRBAAAAChJFAEAAABKEkUAAACAkkQRAAAAoKSOoVzc3t4ejUZjuO7lNDfeeGPKTjXt7e0xduzYlK3nnnsuZaeiQ4cOpf37feqpp1J2qjl48GBs3rw5ZWvZsmUpOxVNnjw5br/99pStJUuWpOxUs3///njyySdTttavX5+yU9GMGTPigQceSNm6/vrrU3aqGTt2bMydOzdla8aMGSk7FTUajbjmmmtStvr6+lJ2qunv748///nPKVtf+9rXUnb+HW+KAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJoggAAABQkigCAAAAlCSKAAAAACWJIgAAAEBJHUO5uK2tLTo6hvSP/J8tXbo0ZSciYt26dWlbrdbT0xMf/ehHU7amTp2asjNo9+7dqXutNHny5Fi8eHHK1h/+8IeUnYiIffv2pW21WltbW4wcOTJl65vf/GbKzqDly5en7rXSjh074u67707ZWrlyZcpO9larTZs2Le2M7Ny5M2Vn0LZt21L3Wmn//v3xy1/+MmXr5ZdfTtmJiLjpppvStlrt+PHjsWvXrpStJUuWpOwMeuihh1L3WunEiRPR19eXsrV3796UnYiIa6+9Nm3rv0FbW1vKzqOPPpqyExHxjW9844w/86YIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJQkigAAAAAliSIAAABASaIIAAAAUJIoAgAAAJTU1mw23/rFbW2vR8T/DN/ttMw7m83m5FbfRIZz+BlGeI7nAs/w3OA5nv08w3OD53j28wzPDZ7j2c8zPPud8RkOKYoAAAAAnCv89RkAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgJFEEAAAAKEkUAQAAAEoSRQAAAICSRBEAAACgpP8FRmTg+2kRNcIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd7RU1RXH8YMJChGkSRGQqqAgIooFaygi0SgEltFgAVyoKMsoaLAbXVl2JWjE2EBYKpa4Is2CFU3UiEKIBQEjRZSHiihgQUlC/sjK9rc3b56Umcd9M9/PX3s4h+Eyd26Zu87eu9qGDRsSAAAAAAAAsmW7bb0BAAAAAAAA2BgPbQAAAAAAADKIhzYAAAAAAAAZxEMbAAAAAACADOKhDQAAAAAAQAbx0AYAAAAAACCDfrw5k6tVq0Z/8G1kw4YN1fLxPuzDbWrlhg0bGubjjdiP2w7HYlHgWCwCHItFgWOxCHAsFgWOxSLAsVgUyj0WWWkDVJ6l23oDAKSUOBaBrOBYBLKBYxHIhnKPRR7aAAAAAAAAZBAPbQAAAAAAADKIhzYAAAAAAAAZxEMbAAAAAACADOKhDQAAAAAAQAbx0AYAAAAAACCDeGgDAAAAAACQQT/e1huA0lGvXj2LO3fubHG/fv3cvObNm1tcs2ZNi6tXr+7mrV271uKHHnrIjU2bNs3idevWbeEWY1v70Y9+VO6f//vf/67kLSlN2233/XP9Fi1aWLzDDju4eUuXLrWY4w0AAGBjNWrUcK+//fZbizds2FDZm4MqhJU2AAAAAAAAGcRDGwAAAAAAgAzioQ0AAAAAAEAGUdMGBbPbbru510OHDrV4wIABFtevX9/Ni7Vr/i/WN9G6Jp06dXJjixYtsnj27NmbuMXYFqpVq2bxiSee6MZuvvlmiydNmmTxJZdc4uZ99913Bdq60qL7IqWUevbsafE111xjcZ06ddy8mTNnWnzddde5scWLF1tMvnbl0FpEKflz57/+9S+L2R+lRWtR/ec//7F4/fr122JzAKBo6HW2bt26bmy//fazuHfv3m6sQ4cOFt91110Wz507181bsmRJPjYTVRgrbQAAAAAAADKIhzYAAAAAAAAZRHoU8qpVq1YWjxs3zo21a9fOYl2iv3z5cjfvjTfesFiX7++///45369x48ZurFevXhaTHpVtmpLTsmVLN6Zt3V9++WWLSYcqjFq1arnXw4YNs3iPPfawOLasbNasmcUxdWrIkCEWf/3113nZTvyPHjt77rmnxUcccYSbN2vWLIt1ybWmmG7Ov6XLwHfccUc3T1NwPvvsMze2Of8etoymxp177rlu7KKLLrJYvxMDBw508/S8i+/9+Mff3zLrvU5KKXXr1s1ivTdZuHChm/f4449bvGrVqjxvoaffhbi9jRo1sljPCTFVjmN2Y3ouJMW0+MQSDXpN0+MmpZT69u1rsV53NR0qpZRq165tcSz1oOcVPWaff/55N++mm26yWNNbUTpYaQMAAAAAAJBBPLQBAAAAAADIoEpJj9KlhLqUukWLFm7eypUrLf7888/dGB0vsil2m9Elwk2aNHFjuuzv/ffft3jUqFFu3jvvvGOxLgGM6VETJkywuF69em6MJb1Vhy4V7dy5sxubN2+exS+++GKlbVMp0WO4T58+bkyPOV0yHI973YcHH3ywGzvqqKMsnjx5ssWcx7eedqg455xzLI7Lr6dMmWJxPs6NulxcU1FTSmnXXXe1WM/RKaX0xRdfbPW/jYppeuJll13mxvQ6ue+++1rMsVi+mCbRvn17i2+99VY31rFjR4v1/vWee+7J+f75TrOJ26vHZuzM+NFHH1m8dOlSi1esWLHV25El8Vq1qZ/z9ttvb3GbNm3cWIMGDSzWzy6mg37zzTebvJ0oPE1DSsmnef/kJz+xeNCgQW6e3sPEtCd9j3Xr1lkc05f0t2/cDp2r909/+MMfyvlfoJSx0gYAAAAAACCDeGgDAAAAAACQQTy0AQAAAAAAyKBKqWmjebbagvKMM85w87TmyZo1a9yY1jmZNm2axa+//rqbp3m6ml8YaQ5hzHnVXFbdplgLQNsOl2oNlfjZaXvtRx55xI3ttttuFo8dO9biOXPmuHlav0jfX/88JZ+7r/sppZTmz5//g9uObNh9990t7tKlixt79dVXLaYeRmFoXZTzzjvPjWldKq2TEusCaI52rGWlx3pZWZnFr732mptHXY0fFs9z/fr1s1hz7eO595NPPtnqf1v3j17vtIZNSr7uR2yPyjFceMcee6zFO+20kxvT66neU3311VeF37AqKNae6NSpk8UNGzZ0Y3p+fO655yyOx6Le2+b7nKf3RCn52hz77LOPG5s5c6bFX375ZcG2aVvbnP+P1jW58MILLR4yZIibp/ta9+fixYvdvJEjR1rMPWnl0d+cbdu2tfjAAw908/T6pNdSrfcV3y/69ttvLZ47d67FDz74oJun18lYO7B+/foWL1myxOJly5a5ecV2bOaDXtP0+D311FPdvJ/97GcW63ciJX+dXLhwocU33HCDm6fnTN3vlYmVNgAAAAAAABnEQxsAAAAAAIAMqpT0KF3S9dZbb1kclxJqu19to5aSX77UtWtXi2Nr01q1alms6UtxKZO2ZGzcuLEb0yWxNWvWtDi2cNN2xCeffLIb05bWxSx+JosWLbL4jjvucGO6jE3bu8e0J6X7d8yYMW5M25euX7/ejb388ssVbTa2oZhSN3ToUItjOsVTTz1lcfyuYcvpcXX55ZdbrC2+U/LLgjUlpqJ9EfevHqd33323xcccc4yb98EHH/zQZpe8mO7Sv39/i/Vznz59uptX0Tl2a+29997utaY7VrSsHPkRU3guueSSnGN6DOuxyLL78sXjRl/Hc+CkSZMsvu666yyOLbTzfR3Tc/ngwYPdmKYbr1692o3p9q5duzav21RVxHTTU045xeIRI0ZYrL8DUvL3m3pObtq0qZunLZvjvtFSDth8uu+0BXtK/t5C0570d0dKPi20du3aFsffE3p8aIv3lFK65ZZbLH788cfLfe+U/Dn2pptucmP6u3XVqlUWxxIfpXqe1ucBe+yxhxsbPny4xdqaPabpx2cFuejfO+CAA9yYpkfdeOONbkzT/fXZQ773GSttAAAAAAAAMoiHNgAAAAAAABnEQxsAAAAAAIAMqpSaNpofqDl/2nIypZTatWtnsbbuSimlgw46yGKtexHz21q3bm2xtu6OtQA0RzXmLWvOsc6LOXFa/+Giiy5yY8OGDbO4lNqB62e5fPlyN6Z5iZrzF2muqrYvjW34tIaDtqxMydcsQrbocZlSSkceeaTFsR6K5pAif371q19ZfNZZZ1kca5BoPq4e27E+wscff2xxzC/XXHFte3nFFVe4ebodMae8lOkx0aFDBzfWvHlzixcsWGBxoWuq7bjjjhbvtddeOcfieRn516pVK/e6RYsWOedqy3W9F0P5Ygvtww8/3OLYvv6FF16wWGtnFLoWm9aQGjVqlBvT+9d7773XjVH3b+O6FxdccIHF+tl9/fXXbp7eX2q9zF122cXN098IWkMqpZROOumkct8P5Yu/v/Ta16tXLzemrb3vv/9+i1988UU3T6+t119/vcXx92JZWZnFsVbNltSLi/dP8XUpivWltJ7Mr3/9a4s7derk5uk9pYrHrJ6v4/2lPm+I+1717NnT4tg+furUqRZrTbNYu3drnwew0gYAAAAAACCDeGgDAAAAAACQQZWSHqV0aVBcwq2vY5rElClTLNalTLoUOyW/NF+X033zzTdunqbnxHaHuiRWl8xp67i4jbENHO2JN251Ftuu/19c9qjL326//XaL4/I5ff8TTjhhi7cTlatu3bru9W677Wbxe++958Y+++yzStmmYle/fn33esyYMRZr2mKk57E333zT4osvvtjN0/Sohg0bujFdct6tWzeLBwwY4Obp+UGXw6ZUWimmkS7X1RSylPy1SpffFzq9rEePHhZ37NjRjWmLY21fivzRe48jjjjCjenxHO9DHn74YYvj8nH8j7ZJP+yww9zYz3/+83LnpbTxda2Q9JzwzDPPWBzPvbo0/84773RjpXpO1fvN2NJX06X0N8OTTz7p5k2fPt1iTf8dOXKkm9eyZUuL43dp6NChFmsL6FJt61weveeP6Sj6eyymLE2YMMFibcW8qd/5WNoBhaG/5S+77DI3piVG9Fwb26DPnz/fYn1O8OCDD7p5FaUgalkV/XePP/54N0/Pu/F8r9/PQw891OL4XdLt35LnBKy0AQAAAAAAyCAe2gAAAAAAAGRQpadHbaq4RFArdK9Zs6bcOCVf5XtL6TKqv/3tbxbr0tiU/BJ0rRydEkscy6OfiS7vjsvMevfubbEun4tLyZYsWWLx888/n6/NRIHp/k3JL32MHS1KdQl3PugxdsMNN7ixmC71f/EYe+mllyzWdKbY7UCP7ZjuqOka++23n8U777yzm3faaadZHFMpzz///HL/rWIUU4OPO+44i2OXDE0nLHQ3GN2vV199tcWx06N2jIppycgP7cB35plnujH9/sSl5JrujfJpyr2ek1JKqWnTphbHz1Y7mWrKTEwR3JLzVywDoClR2kUnXi+1SyAdiv5Hz1ennHKKG9Przpw5cyy+8cYb3TxNO9P07vj56zkz7sM+ffpYPHr06JzvUcoaN25s8a233urG9HfDtdde68Zmz55tMZ9ndsRU/Mcee8zi7t27uzE9T2rnp0ceecTNGzt2rMX6mzBXOY6UNr5H1TIM9erVs7hWrVpunv5WiffKes7XznPxfm5rS6ew0gYAAAAAACCDeGgDAAAAAACQQTy0AQAAAAAAyKDM1rTZljTv7phjjsk57y9/+YvFCxYsKOg2FYNcdS923XVXN69NmzYWa15irF+kOcG0WM82zQWNbYu1XtUTTzzhxoq9fkkhac73L37xCzemebb6GWt7zJT8Mfbdd99t0r+r+zOllF5//XWL//rXv1rct29fN2+HHXaweNCgQW5s/PjxFr/99tubtB1VlX4OKfnPItaP0Va0scZGvrVt29ZibZEZjRs3zmLqCRSG1qTSWirRypUr3Wta2f4wrcMUjze9b9G6BSn5Wm3vvvuuxXr+S8nXSdD9E8+v2i76j3/8oxvr0KGDxXr+njZtmps3a9asBE/3W/Xq1d2Y7g+tU7lw4UI3T69xsd10LrG2hZ4btbU158zv6fdc45RS+vrrry2O57VNvVdB5dL6TymldPDBB1usvxFS8sfYvHnzLH7ggQfcvGXLllmsx058P32ttRVT8vdYRx99tMVaOy6Kx6nWuZo8ebLF+j3NB1baAAAAAAAAZBAPbQAAAAAAADKI9KjklyamlNIZZ5xhcZcuXSxesWKFm3f22WdbXFF7MWxMl5399Kc/dWOdO3e2WJcqT5kyxc3T5WjItiZNmljcqVMnN6ZLkl955ZVK26Zi16xZM4vjUn5dUq9to2M79nwsM9aWt9ratEePHm5enTp1LI6tFvWcUOzpUfH/rqkwn376qRt7+OGHK2WbUkrp0ksvtVjTRNauXevmaUtO0hsL48ADD7Q4thLWz1xT1VLaOHURG1u/fr3F06dPd2N6HorL73U/jBw50uKY/q1/T6998TjaeeedLW7YsKEb01QbvQ8aPHiwm8fxtzFNV1i6dKkb032obYDjMaapqJrOqtfSlFLaZZddLI7H3ocffmhx/A1SyvSzuP766y2O9zB6DYrXTNLNqgY910Z6jmvevLnFw4cPd/PmzJljsR6XPXv2dPP22msvixs1auTG9LsV09OVfpfi788RI0ZY/Mknn1ic73MwZwoAAAAAAIAM4qENAAAAAABABpEelTauSj5q1CiLdanUDTfc4OZ98MEHhd2wIqNLFi+88EKLzzvvPDdPK/rr8uSLLrrIzcvisse4zLWUu1rp8sbu3btbHDtyaErUF198UfgNK1KxO4We1+JSfk17uvLKKy3+8ssv875dujx00aJFFmvaVEp+iXNcSq5L1XN1vioWcan3559/bnFZWZkb088l3+Jxuu+++1qsHVOee+45Ny92LEJ+6LVlyJAh5f55Sv7Yvvvuuwu/YUVGzyl33HGHG3vhhRcs1iX7KaV06KGHWqxdR2NXN035fvzxxy2O95OaqqrvnZLvpNm/f3+LV69enVAxPT7+8Y9/uDFNodDyB4cddpibp/tQz8/x3K37Kaa/aUcqTfUpdXr8VZSqovc0V111lRt75513LH7//fctjunexXj/kGUxfVDPr0OHDnVjeixpeYXjjjvOzevXr5/Fun/j/YuW5Ij3yrnEsif6W2XYsGFu7J///KfFhfzdx0obAAAAAACADOKhDQAAAAAAQAbx0AYAAAAAACCDSramjeZDalu5lHyrRc0zvu+++9w88iE3j7a+1JaYMffwo48+sviss86yWFs1ZlUp17CJNB/53HPPtTgeN7fffrvFWaxTVFXEPN1jjz0251xtHf38888XbJtS8ufajh07Whzz1XXfa82UlHx71GIX61JojZjY+ldbymrNhE09D8XvjOZ9H3744W6sVatW5b7HQw895F5zDBdGgwYNLD7yyCNzzluxYoXF2noUm09rl6SU0ty5c8uNU/L19y699FKL4zGW6/ioUaOGe621inbffXc39sQTT1is9Tvww7Re2ltvveXGtBVw69atLT7wwAPdPK1TpMdY/L4sX77c4rjfte7R008/bfHf//53N6/Ufmfo/1fbOcfao3p8tG3b1o09/PDDFr/55psWv/32226e3vto+/dYy0T3azxO9T5Gvxfx90pF7a1LRfxcL7/8cosnTZrkxnr06GGxHn+NGzd28/T8qr8l27Vr5+bpfWisA6ffOd3XsSac1lWN9ZEqCyttAAAAAAAAMoiHNgAAAAAAABlUsulR2trvkEMOcWO6fPKcc86xuBDtcItZTH/Q5W+6jC0uG9TPPLYFRtWhy/nbt29vsbbITCmll19+udK2qZjFtqHaojku0dfW6nFJ99aK7cX33ntvizUtsm7dum6ebqMuVU7JLzMv9uXisTWsLuk+/vjj3djFF19s8ezZsy1u0aKFm1evXj2LNcUqfpbahvaggw5yY9qCU5d+azpOee+J/DjzzDMtrlmzpsXx89bUAFLVto0t+dxjSqO2+Y7porfddpvFer+KH6afs6bOpJTSvffea7GmQtSuXdvN02uVXkufffZZN09TKPQ6mFJKBxxwgMVaouH000938/RaWGrn1t/97ncW6/1MSj5dN6Ys7bnnnhZrWlU8xjQ9R1ObysrK3Dy9R4pt3dWSJUssnjBhghubNm1aue8Xt0u/W8W+v/U8GVPX9LV+JvH+Un9n9uzZ0+JY9mSnnXayOP7m1HuuUaNGWTxx4kQ3LwvlL1hpAwAAAAAAkEE8tAEAAAAAAMggHtoAAAAAAABkUMnUtKlevbp7rW2GYy2IP/3pTxZrKz5snuOOO8693meffSzWXMannnrKzdO84CzmdMb6IBXloGZx+ytL9+7dLdYaRosWLXLzYo0bbJmY1605vPE7q2OaK/7qq6+6eZr7q9/l+H477rijxSeccIIbO+200yzWc0CseaX5/zEfXOsGFLtYo+LBBx+0uFu3bm5M913//v0t1tbd8bWeez/++GM3T/dxrDmkY7qNcT+WUk5+IcXcfa1po2L9lLvuustiPv+qI7Zx15onb7zxhht7//33K2Wbil2sLTJ69GiL9T701FNPdfP0HPfaa69ZPGvWLDdPf3fovJT8747OnTtb/Mgjj7h5+m8vWLDAjRX78a3f83hfcdJJJ1k8cOBAN9a0aVOL9doX71t0H+j9k9aAS8l/zvEz12ths2bNLG7ZsqWbp//29OnT3Zi2wq6oRlWx7+9c9P8d69Ho9a9r164W675IyX/+a9ascWO/+c1vLL7vvvsszkINm4iVNgAAAAAAABnEQxsAAAAAAIAMKur0KF0OdfLJJ7uxjh07WhyXmp5//vkW0zJz82y33ffPAQcMGJBzni7Lv/HGG92YLhUs9FJ7ff+4dFL/L5oqEFMU1MyZM91r/f5om9xipJ9XSikNGjTIYv1sp06d6uZxjOVH/BwrWhac6/scz4XaClHfr127dm7eiBEjLNa2iyn5FpnxO6I+/PBDix944AE3lsVlqpVF98mQIUPc2P7772/x4YcfbnFcmq3pibpM//7773fz9Nx72WWXubH27dtbrOeylStXunmluoQ732J6WqNGjcqdpy2By3uN7GrSpInFMf1Nz7fz5893Y7T5Lgy9zsydO9fimJak6cV6zoz3ePp+ixcvdmMfffSRxZqSsccee7h5ei08+uij3VhMby02ei2ZN2+eG/vtb39r8c033+zG9H7kvPPOsziWwtD7kVhCI5eK0nP0mG3cuLGb17dvX4ufe+45N/bVV19t0r+Nje29994W62/3mLatx+Ytt9zixvQYy/q9JittAAAAAAAAMoiHNgAAAAAAABlU1OlRDRs2tPiqq65yY7qk7brrrnNjcbk3Np0uDzzssMPcWK7UiDZt2rjX7733nsXaUUY71KTkK/+vW7fOjekSN90mrSqfUkqtWrWyWLuvpOTTOrp06WJxTDXR5ej6d1LynV+KXfy/77XXXhZ/+eWXFk+cOLHStqmUxGW72nFp5513dmM1a9a0uHfv3hbHc5929jriiCMsjp0c6tevb3FFKVC63HnVqlVu7MQTT7R49erVOd+j1OhnFpfDP/HEE+XGFb1HRelLunz8k08+yfkeel4uKyvL+X7YcjHNULtJ6f3L1Vdf7ebF8wCyRe8fOnXqZLGm7Kfkj7d4TtXXpBcXnqZApeSvXVt6vL377rsWa3rUJZdc4ua1bt3a4t13392N6Tm61NJSNUUwdiDV33R6fMQURO0YpfeoMbVGf0/ovJT8NVN/o8TfCfo9iemNufZdqe3TTRG7Ks6YMcNiva+N58Unn3zS4pgeVZXSTVlpAwAAAAAAkEE8tAEAAAAAAMggHtoAAAAAAABkUNHVtNF8t7PPPtvievXquXnaQrGiWgDYPLVr17ZY80VT8rmfWmPj4osvdvOGDx9usbau1ba1Kfk2edouOKWUPv30U4tbtGhhsdbeSMl/L7T2TUq5c8W1nkNKPo81tjkvpZo2PXr0cK91H2vuNTUwCiPm5Y4bN87iK6+80o1pzvYhhxxisdZuSsnncuvxHHO+K6pjo9u1aNEii0866SQ3b86cORaTy71pCtmeMl4zlbar1bbwyJ9YY01rJOhnPnXq1ErbJhROrJWh9Rk0Tsm3nNb6Kpw3CyOeZ/Nx3tX7SK1pE+8hGzVqZHH79u3d2KxZs8p9v1ITv/f62+D3v/99znlaN0zrW8bfAtouOtbP1JqWWp8zzhs/frzFsWYfx23F9Nr32GOPuTGtXauf42effebmjRw50uJ4rq1KWGkDAAAAAACQQTy0AQAAAAAAyKAqnx4V26odfPDBFp922mkWx+VnuhyR5d2FsWLFCvdal/Tq8sO2bdu6ebpPK0q70LSkBg0a5HwPFVOedJlr/B7oa93e+B533323xboUsxTo53zFFVe4MU1V1O+CtmpH/sRznKZHxRbdnTt3trh69eoW6zEa31P3da7jK6WNl2m/8cYbFg8cONBiXVaMbMjVvjQln+am58Cq1C4z6/S4iu199TNfvHixxXGpPbJN72kaN25ssR57KfnzaLzn0JSANWvWWFxRK+GKztmkZ1Q+/cw1xU3TaFJK6fTTT7c4pkdpy/i5c+eW+97xdantaz0/3nbbbW5MP4vBgwfnfA+9FsZjTEtCzJs3z2K9/0oppddee83iXKUXUL5f/vKXFh999NFuTM9rum/GjBnj5mlKd1XGShsAAAAAAIAM4qENAAAAAABABvHQBgAAAAAAIIOqfE2bZs2audejR4+2uGnTphYvXLjQzbvnnnssLrUcz0LSVtu9evVyY3369LF4//33t7h169ZuntbV0DbfWiMlJZ/z/d5777kx/XtvvfWWxbFV40svvWTx559/7sa0LVy3bt0sju3FX3311XL/TinQXF9tmZiS/6wnT55c7p+jcLTl4VFHHeXGHn30UYsPOOAAi2Ory1w1EWJO9sqVKy0eNWqUG/vzn/9sMfWMsk1buXft2tWNae2j2IIY+aHHWLzeab6+Xu84n1Ytul/1mlm3bl03T+9h9tlnHze2yy67WFxWVmZxPL+uX79+q7YVlUPvZRctWuTG9Frbt29fN6atqLV2S6wnqfNK+feO1g5KKaWpU6da3L9/f4tjjUyln2VKKT399NMWz5gxw+I333zTzePeZ/PoPcbYsWMtrqjOqdatufXWW91YsXzvWWkDAAAAAACQQTy0AQAAAAAAyKAqmR6ly0uHDh3qxtq1a2fxt99+a/Gll17q5sUlbsg/XbabUkr33ntvuXFVsGDBAouLZZldPuiyXk2DSckv946tFlG5NH0pJZ+62KRJE4u7dOni5un5VM+ZzzzzjJu3bNkyi/W8i6pFU6BiuoYuS6ZlaWFo2+e4nH7dunUWx/QHVB2a5lavXj2La9SokfPvxGOxQ4cOFut91vLly908/T7FVCmO4ezQe0q910zJ79+Ylqqt3/X7o+eKlEih/L94765t0o899liLe/To4ebpcTpr1iw3Nn/+fIs1RS1+5vxu2DyaClinTp2c8/Qz79evn8VfffVVYTZsG2OlDQAAAAAAQAbx0AYAAAAAACCDqmR6VPPmzS0ePHiwG9Pl3dph4dlnny34dqF4sbSxfJoeFVMVkV26RF+7ocXOaCgta9assfiuu+5yY4MGDbL4wgsvtJh0uPzRFJZrrrnGjWl3E+2MwbWpatG0JN2PMQVKO0bFjkLapfLjjz/O+W/VqlXLYk0jSMlfA5AdMZVZuxDFlLaZM2du0ntqV8D4HnwP/kdTvCdOnLgNtwQp+eudpnnGtLPx48dbHDt2FSNW2gAAAAAAAGQQD20AAAAAAAAyiIc2AAAAAAAAGVRlatpUq1bN4kMPPdTiBg0auHma362t2WL7TAAA8D2tb3Dttde6sfgahTVlypQKX6PqW7p0qcWnn366G9M6JLFdN3VIin4T6RgAAAHESURBVFesOXPnnXdaHOsGap0i/S7F7wffF2Sd1q1Jybe31/bdr7zyipt3+eWXW1wKre1ZaQMAAAAAAJBBPLQBAAAAAADIoCqTHpWLtiiNrx999FGLaYsJAACArIlL+0npL03fffede11WVmbxLbfckvPvrV27tmDbBBTa9ttv715/8cUXFs+YMcPiCy64wM3T1KlSwEobAAAAAACADOKhDQAAAAAAQAbx0AYAAAAAACCDqm1OrZdq1aplrjBMnTp13GttE7ZixYrK3pyC2bBhQ7UfnvXDsrgPS8jsDRs2dM3HG7Eftx2OxaLAsVgEOBaLAsdiEeBYLAoci0Wg2I7FatW+/++UUH3aco9FVtoAAAAAAABkEA9tAAAAAAAAMmhzW36vTCktLcSGbKnVq1dX+LpItMzje2VuH5YQ9mPVxz4sDuzHqo99WBzYj1Uf+7A4sB+rvqLbhyWUEqXK3Y+bVdMGAAAAAAAAlYP0KAAAAAAAgAzioQ0AAAAAAEAG8dAGAAAAAAAgg3hoAwAAAAAAkEE8tAEAAAAAAMggHtoAAAAAAABkEA9tAAAAAAAAMoiHNgAAAAAAABnEQxsAAAAAAIAM+i+9tKNOsGxoCAAAAABJRU5ErkJggg==\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: 2.66 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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0034, PSNR (dB) = 24.63\n",
      "MSE = 0.0057, PSNR (dB) = 22.47\n",
      "MSE = 0.0009, PSNR (dB) = 30.59\n",
      "MSE = 0.0068, PSNR (dB) = 21.64\n",
      "MSE = 0.0049, PSNR (dB) = 23.08\n",
      "MSE = 0.0008, PSNR (dB) = 31.04\n",
      "MSE = 0.0091, PSNR (dB) = 20.40\n",
      "MSE = 0.0108, PSNR (dB) = 19.67\n",
      "MSE = 0.0132, PSNR (dB) = 18.80\n",
      "MSE = 0.0061, PSNR (dB) = 22.16\n",
      "Wall time: 10 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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9624\n",
      "0.9449\n",
      "0.9864\n",
      "0.9403\n",
      "0.9312\n",
      "0.9907\n",
      "0.8899\n",
      "0.8784\n",
      "0.8349\n",
      "0.9470\n",
      "Wall time: 274 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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9921\n",
      "0.9832\n",
      "0.9980\n",
      "0.9851\n",
      "0.9846\n",
      "0.9982\n",
      "0.9759\n",
      "0.9707\n",
      "0.9621\n",
      "0.9828\n",
      "Wall time: 898 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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Deep Autoencoders in Keras - 784-256-128-64-32-64-128-256-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
}