Deep Autoencoders in Keras - 784-256-128-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 = 64"
   ]
  },
  {
   "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) # the bottleneck = encoding_dim\n",
    "\n",
    "decoded = Dense(128, activation='relu')(encoded)\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[-3](encoded_input)\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": [
      "1595257900\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 51s - loss: 0.6908WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.212496). Check your callbacks.\n",
      "235/235 [==============================] - 84s 356ms/step - loss: 0.2134 - val_loss: 0.1369\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 63s 269ms/step - loss: 0.1230 - val_loss: 0.1110\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 66s 280ms/step - loss: 0.1072 - val_loss: 0.1008\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 66s 280ms/step - loss: 0.0992 - val_loss: 0.0955\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 79s 334ms/step - loss: 0.0945 - val_loss: 0.0916\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 75s 321ms/step - loss: 0.0913 - val_loss: 0.0896\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 78s 332ms/step - loss: 0.0890 - val_loss: 0.0871\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 88s 374ms/step - loss: 0.0870 - val_loss: 0.0861\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 84s 358ms/step - loss: 0.0855 - val_loss: 0.0842\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 89s 380ms/step - loss: 0.0842 - val_loss: 0.0831\n",
      "Epoch 11/50\n",
      "235/235 [==============================] - 91s 386ms/step - loss: 0.0832 - val_loss: 0.0833\n",
      "Epoch 12/50\n",
      "235/235 [==============================] - 91s 388ms/step - loss: 0.0822 - val_loss: 0.0811\n",
      "Epoch 13/50\n",
      "235/235 [==============================] - 88s 374ms/step - loss: 0.0814 - val_loss: 0.0811\n",
      "Epoch 14/50\n",
      "235/235 [==============================] - 90s 384ms/step - loss: 0.0808 - val_loss: 0.0799\n",
      "Epoch 15/50\n",
      "235/235 [==============================] - 75s 321ms/step - loss: 0.0801 - val_loss: 0.0791\n",
      "Epoch 16/50\n",
      "235/235 [==============================] - 76s 322ms/step - loss: 0.0795 - val_loss: 0.0792\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 66s 282ms/step - loss: 0.0790 - val_loss: 0.0781\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 65s 277ms/step - loss: 0.0785 - val_loss: 0.0778\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 64s 274ms/step - loss: 0.0780 - val_loss: 0.0777\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 65s 277ms/step - loss: 0.0777 - val_loss: 0.0772\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 60s 254ms/step - loss: 0.0773 - val_loss: 0.0768\n",
      "Epoch 22/50\n",
      "235/235 [==============================] - 55s 233ms/step - loss: 0.0769 - val_loss: 0.0766\n",
      "Epoch 23/50\n",
      "235/235 [==============================] - 48s 205ms/step - loss: 0.0766 - val_loss: 0.0762\n",
      "Epoch 24/50\n",
      "235/235 [==============================] - 56s 236ms/step - loss: 0.0763 - val_loss: 0.0758\n",
      "Epoch 25/50\n",
      "235/235 [==============================] - 59s 250ms/step - loss: 0.0760 - val_loss: 0.0756\n",
      "Epoch 26/50\n",
      "235/235 [==============================] - 59s 253ms/step - loss: 0.0758 - val_loss: 0.0755\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 53s 225ms/step - loss: 0.0755 - val_loss: 0.0749\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 70s 299ms/step - loss: 0.0753 - val_loss: 0.0750\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 50s 215ms/step - loss: 0.0751 - val_loss: 0.0747\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 45s 191ms/step - loss: 0.0749 - val_loss: 0.0747\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 65s 275ms/step - loss: 0.0747 - val_loss: 0.0746\n",
      "Epoch 32/50\n",
      "235/235 [==============================] - 50s 211ms/step - loss: 0.0744 - val_loss: 0.0746\n",
      "Epoch 33/50\n",
      "235/235 [==============================] - 54s 232ms/step - loss: 0.0743 - val_loss: 0.0740\n",
      "Epoch 34/50\n",
      "235/235 [==============================] - 43s 185ms/step - loss: 0.0741 - val_loss: 0.0740\n",
      "Epoch 35/50\n",
      "235/235 [==============================] - 54s 228ms/step - loss: 0.0740 - val_loss: 0.0738\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 56s 238ms/step - loss: 0.0739 - val_loss: 0.0738\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 59s 250ms/step - loss: 0.0736 - val_loss: 0.0736\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 53s 224ms/step - loss: 0.0735 - val_loss: 0.0733\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 58s 245ms/step - loss: 0.0734 - val_loss: 0.0732\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 58s 247ms/step - loss: 0.0733 - val_loss: 0.0731\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 47s 199ms/step - loss: 0.0731 - val_loss: 0.0727\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 55s 234ms/step - loss: 0.0730 - val_loss: 0.0731\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 49s 208ms/step - loss: 0.0729 - val_loss: 0.0728\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 54s 230ms/step - loss: 0.0728 - val_loss: 0.0727\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 58s 245ms/step - loss: 0.0726 - val_loss: 0.0726\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 56s 236ms/step - loss: 0.0726 - val_loss: 0.0726\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 39s 167ms/step - loss: 0.0725 - val_loss: 0.0725\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 39s 164ms/step - loss: 0.0724 - val_loss: 0.0723\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 38s 162ms/step - loss: 0.0724 - val_loss: 0.0723\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 45s 193ms/step - loss: 0.0722 - val_loss: 0.0722\n",
      "1595257900\n",
      "Wall time: 52min 22s\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": 8,
   "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": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAATsUlEQVR4nO3b+3/P9f/H8ef7Y5tttGFL04aZyy5ConFREsVFUk455NDJqQuSuCRUaA5LISOHhIt0ySFUhDnORTHHRDmGObTPGDF2CGMOe33/gC/reb9c3n32rMvt+mOX2x57ba/36/V+vR+Wz/M8AwAAAAAAALf8p6QPAAAAAAAAAP8fSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAcFKHFoaKhXrlw56z4oKEg7mADpcMzly5et24oVK0qzc3JypL5q1apSf/r0aev22rVrprCw0Cd9g3uIjIz0YmNjrfvs7Gxpfm5urtSXLl3aus3Pz5dmR0ZGSn1hYaHUK8dujDHnz5+/5Hne/dIX3UOFChW8mJgY61792a5cuSL1ZcqUsW59Pu2lrJ73qKgoqVeu9ZycHHPt2jW/XIvqOTx37pw0PzQ0VOrz8vKsW+V9wBj9vlC2bFmpV1/fubm5frsWw8PDPeU1d/36dWm++r74n//Y/1vMzZs3pdl37tyR+rCwMKkPCQmxbjMzM82lS5dK5FosVaqUNP/ixYtSr7z+lfNtjH4/LV++vNQr9xFjjDl37pzfrsXQ0FAvPDzcug8ODpbmq/cxz/OsW/Weqt4X1POYlZVl3ebn55uCggK/XIshISGect9QX/+qiIgI6/bMmTPSbOW5yRhjbty4IfVxcXFSv2/fPr9di+XKlfMqVapk3av3DfW9SLkW/+7nFeV9zhjt/SY3N9dvz6jBwcGe8rOpz/ZFRUVSf/v2betWee0Zo3/WVd871J3CjRs37notSnf9cuXKmT59+lj31apVU8bLbyqLFy+2bgcOHCjNXrp0qdTPnj1b6rt3727dbty4UZpdnNjYWLN3717r/vPPP5fmr1ixQj4eW+vXr5dmv/7661KvLNKM0Y7dGGOSkpL+K31BMWJiYszatWut+xMnTkjzt27dKvUNGjSwbtVlrnrehw8fLvXLli2zbqdMmSLNLk5MTIxZt26ddZ+YmCjNr1u3rtQrr6c2bdpIs9X7wpNPPin1J0+elPply5b57VqMioqS7pNHjhyR5qvLZ2VZp37IUB9mW7ZsKfW1a9e2bps2bSrNLo56LaoP6dOnT5d65WdT//EgNTVV6jt27Cj1KSkpUj9y5Ei/XYvh4eGmZ8+e1n2NGjWk+cuXL5d65UNGu3btpNkVKlSQ+s6dO0v9Bx98YN1++eWX0uzihIWFmW7duln36j9OKB/ejTGmd+/e1u2gQYOk2Q0bNpT69PR0qVc/x/h8Pr9di5UqVTLz58+37tX7hrrkURYEyjOhMcY89dRTUl+rVi2pVxa66ntNccqWLWvatm1r3auL5GvXrkm9cs5HjBghzZ47d67UP/TQQ1K/cOFCqT927Nhdr0X+9ygAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAABwUocWBgoHnwwQet+8OHD0sH06xZM6kvKir622ZPmjRJ6n0+n9S3atXKur1z5440uzhnz541Q4cOte5DQkKk+XXq1JH6n376ybrNysqSZqu/t4ULF0r9hAkTpN6fLly4YJKTk637bdu2SfMbNmwo9Zs2bbJuZ8yYIc0+dOiQ1EdHR0t9UlKSdavcc/5KZmamGThwoHWfkJAgzT9w4IDUt2nTxrodPHiwNLtXr15SP378eKkfOXKk1PvTlStXzI8//mjdR0ZGSvMPHjwo9UFBQdbtuHHjpNmvvPKK1CcmJkr9zp07rdvr169Ls4uTlZVl3nvvPev+5s2b0vwqVapIfXp6unXbr18/afbRo0elvmbNmlJfku+LpUqVMmFhYda98jxrjDGjR4+W+rS0NOt2wIAB0uy3335b6rt16yb158+ft27Xrl0rzS5OcHCwqV27tnWvXCvGGBMRESH133//vXW7ceNGabb6fL1y5Uqpnzp1qtT707lz56Rnq9DQUGl+cHCw1Cvvi3l5edJstf/444+lXnkf9eczamBgoImKirLuK1euLM1Xz3lsbKx126RJE2l2SkqK1Ldt21bqR40aJfUffvjhXf87f2kDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgn+d59rHPZx8bY3r27CkdTH5+vtTHx8dbt7du3ZJm5+XlSf3q1aul/vLly1LveZ5P+oJ7UM/h7t27pflr1qyR+oyMDOv2kUcekWY3adJE6tu1ayf12dnZUm+M2ed5XgP1i+5GPY9btmyR5q9fv17q9+3bZ93Wq1dPmv3AAw9I/ZIlS6T+l19+kfqSuhaTkpKk+QMHDpT63r17W7fqtRgRESH16mukadOmUm9K8Fpct26dNH/v3r3q8Vi3WVlZ0uy6detK/e3bt6X+rbfekvqSuhZTUlKk+WlpaVJ/+vRp6zYsLEya/eijj0r9smXLpH7Hjh1Sb0rwWlSfb8aMGSP16n1Mcd9990n9wYMHpV497yV1LS5cuFCav2HDBqlPSEiwbtPT06XZUVFRUq88ZxmjP4+bErwWe/ToIc3/888/pT42Nta6zczMlGarn+fGjx8v9Y0bN5b6kroWe/XqJc1X74+FhYXWrfpcX1BQIPXqs5D6/G7ucS3ylzYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA4KUOJq1aqZcePGWfcnT56UDqZWrVpS/8cff1i3NWrUkGYHBQVJfU5OjtSvWrVK6v0lKirK9OzZ07rPy8uT5p86dUrqo6OjrdsvvvhCmr1r1y6pHzNmjNS/+eabUu9PZcuWNfXq1bPud+zYIc3PysqS+pkzZ1q3GRkZ0uwTJ05IvXotlpT69eubvXv3WveHDx+W5l+6dEnqb968ad0+/PDD0uyoqCip/yepUaOGmTdvnnU/a9YsaX5ERITUX7161boNCJAeAcxzzz0n9e3bt5f6khIVFWX69Olj3Su/Y2OMKSwslPrvvvvOut2+fbs0+8aNG1Kv3kdKUnx8vJk+fbp1/9FHH0nz1efCxo0bW7cLFiyQZo8cOVLqz549K/UlpXz58qZFixbWfYMGDaT5EydOlPqnn37auo2Li5Nmq5+RunfvLvVr1qyRen+qXr26SU5Otu49z5Pmp6amSn1CQoJ1q17n6j11xowZUl9SoqOjzaBBg6z7wMBAaX5aWprU16xZ07odNmyYNFt5fRij36/9hb+0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAH+TzPs499PvvYGNO2bVvpYFJSUqT+66+/tm4rV64szR41apTUx8fHS/3t27et25SUFHPp0iWf9A3uQT2He/bskeY3bNhQ6hMSEqzbZs2aSbN37Ngh9Y0bN5b6NWvWSP3x48f3eZ7XQPqie1DP49ChQ6X5kydPlvr333/fug0LC5NmBwcHS31qaqrUd+rUybpNSkoyGRkZJXItrl27VprfunVrqX/xxRet26ZNm0qzlfudMcbk5eVJ/datW6V+y5YtJXYtTpkyRZo/ZMgQqd+5c6d1W1hYKM1esmSJ1J85c0bqi4qKrNtdu3aZ/Pz8ErkWe/ToIc3/6quvpH7SpEnWbUhIiDS7oKBA6tVzqN4bunTpUmLXYr9+/aT5c+bMkfpNmzZZt7t27ZJmZ2ZmSr3yrGWMMWlpadbtxo0bTU5OTolci+o5Uc/5tm3brFv12WbmzJlSf+HCBalXj2fRokUldi1+8skn0vxhw4ZJfWJionV769YtaXZGRobUHzp0SOrLlStn3e7fv99cuXKlRK7FVq1aSfM3bNgg9Zs3b7Zu1c+u6vN19erVpf7o0aNSv2fPnrtei/ylDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAgwKUuGbNmmbx4sXWfUJCgnQwXbt2lfpff/3Vul24cKE0+7HHHpP6ypUrS314eLh1u2XLFmn2X33fJk2aWPcNGzaU5v/www9Sr/ye58+fL83u3r271M+bN0/qY2JipP748eNSX5wqVaqYESNGWPf9+/eX5r/22mtS3759e+t28+bN0uzMzEyp3717t9RfuXLFur18+bI0uzhxcXHm448/tu5bt24tzVf7pKQk67Zfv37S7BYtWkh9s2bNpP78+fNS7897av369c3evXute5/PJ81fuXKl1K9du9a6HTZsmDT7m2++kfpOnTpJfZ8+fazbBg0aSLOLExUVJX3v8ePHS/PnzJkj9fHx8dbt9u3bpdk5OTlSv2LFCqlPS0uTen+qWrWqGTlypHXft29faf61a9ekXnkfDQoKkmY3atRI6qOjo6V+yZIl1q0/r8XY2FgzevRo675Xr17S/GnTpkn9sWPHrNvr169LsxctWiT1ycnJUn/gwAGp96eQkBDpPqa+F6n31D179li3L7/8sjT7559/lvqioiKpb9q0qXWrPi8XJzIy0nTs2NG6nzt3rjS/Q4cOUj9x4kTrNi4uTppdqlQpqVd/z+oz7b1er/ylDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAgwKUuKCgwOzfv9+6nzt3rnQwffv2lfpx48ZZty1btpRmx8TESP0TTzwh9crPmpubK80ujs/nM4GBgdb9xIkTpfnNmzeX+p49e1q3BQUF0mzl5zTGmLCwMKkvX7681PtTdna2mT17tnXveZ403+fzSX2dOnWs24MHD0qzO3ToIPXVqlWT+j179li3v/32mzS7OOHh4eb555+37tu3by/NX7VqldQHBwdbt/369ZNmr1ixQurV33P9+vWl3p+ys7Ol97ratWtL81944QWpnzZtmnW7dOlSaXbXrl2lvn///lJ/9uxZ6/bcuXPS7OJcv35dera5ePGiNL9ixYpSP2HCBOt269at0uzBgwdLfXx8vNSnpqZKvfp+UBzP80xRUZF1P3nyZGl+mTJlpH7jxo3WbWJiojRbfeZUr5chQ4ZYt2fOnJFmF+fChQtm6tSp1v348eOl+errX3lWf/XVV6XZ7777rtRXr15d6u+//36pnzNnjtQXJzIy0vTu3du6//3336X56jOI8oywadMmaXZycrLUjx07VuqVzzJXr16VZhcnLCzMPPPMM9a9+vmvW7duUq88O7300kvS7JCQEKmPjY2V+gULFkj9vfCXNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADgpQ4vz8fLN69Wrrvl69evIBKebOnWvdvvHGG9LsFStWSP2hQ4ekfuDAgdZtenq6NLs4nueZoqIi6/7UqVN++953s3//fuv2nXfekWbfvHlT6rt06SL106ZNk3p/CggIMOXLl7fu+/bt+zcejTFBQUHWbd26daXZjRo1kvrJkydL/ZQpU6xb9TovTlZWlklMTLTu4+Li/Pa976agoMC6Ve8Lzz77rNSr19aRI0ek3p9Kly5tqlatat2PHTtWmt+5c2epT01NtW5btGghzVavrRs3bkh9YWGhdet5njS7OD6fzwQGBlr3n332md++991ER0dbt927d5dmq/ffXbt2Sf23334r9cp7x1+5c+eOyc/Pt+4HDBggzR86dKjUDx8+3Lpdvny5NHvevHlSX7p0aakvKWXKlDGPP/64db99+/a/8WiMOXHihHWrPFsbox/7rVu3pP727dtS70+FhYUmIyPDuv/000+l+eozgvL6b968uTR7x44dUq++pyvP+gEB0sf6YmVnZ5tZs2ZZ961bt/bb976b2rVrW7fqtRgaGir1u3fvlvp27dpJ/b59++763/lLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB/k8z7OPfb5sY8x//77DwT1U9Tzvfn8M4hyWKM7jPx/n8N+B8/jPxzn8d+A8/vNxDv8dOI//fJzDf4e7nkdpaQMAAAAAAID/Df73KAAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAH/R/jOnpecEOm3QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de9SVY/7H8W+GSeispKLQgZFJCitjFA2SUxYih5wPoxkxy6lQzguDQRZlJKoRyqmxlsHQYHLOoQOrnIpqOikVOc48vz+s39fn+nr20/M87f10773fr78+e66r/Wz73te9733P9b2uehUVFQYAAAAAAIBs2WhDvwAAAAAAAAD8HDdtAAAAAAAAMoibNgAAAAAAABnETRsAAAAAAIAM4qYNAAAAAABABnHTBgAAAAAAIIM2rknnevXqsT/4BlJRUVEvH8/DMdyglldUVLTIxxNxHDccxmJJYCyWAMZiSWAslgDGYklgLJYAxmJJqHQsMtMGqDvzN/QLAGBmjEUgKxiLQDYwFoFsqHQsctMGAAAAAAAgg7hpAwAAAAAAkEHctAEAAAAAAMggbtoAAAAAAABkEDdtAAAAAAAAMoibNgAAAAAAABnETRsAAAAAAIAM2nhDvwCUj4YNG3ru3r2758svvzzp16ZNG89bbrml50WLFiX9fvnLX3p+4403kraTTjrJ8//+979avmIA/2+TTTbJ2fbDDz94rqioqIuXAwAAkHn16tXz3KpVq6Rt2bJlnvVaCoiYaQMAAAAAAJBB3LQBAAAAAADIIG7aAAAAAAAAZBBr2qBg2rdvnzweOXKk5969e3tu0KBBzufQOtDmzZsnbbp2RqNGjZI2Xe/mm2++qdbrxYahx/hXv/pV0nbjjTd6njx5suf77rsv6cc6KoWx2Wabeb7qqqs8H3fccUm/xYsXex40aFDSNmvWrAK9OlSXjjHGSvn6xS9+4VnXeuMzAQD5s+mmmyaPO3fu7Dn+5tH1Av/73/96XrlyZdJv7dq1+XyJKELMtAEAAAAAAMggbtoAAAAAAABkEOVRyCvd1nvMmDFJ2z777FPpv1m1alXy+JFHHvGsUwX79euX9GvdunWlf9fM7NBDD/U8adKkdb1sZMS+++6bPN55550933LLLZ6Zzl8YG22U3sfv06eP56OPPtpzy5Ytk346/k488cSkbejQoZ61JAP5tfnmm3vu2bNn0qYlakuWLPGcj3EUPzP6WYjTufX467kdhXHYYYclj2+++eZK86hRo+rsNZWK+vXrJ4+1LLRp06aeX3vttaTfiy++6Lkuv8diuUaPHj08z5kzx7NuP4zKUWZY2ho3bpw83n333T3H3zHdunXzrGVPbdu2Tfq1adPGc9zWe+ONf/opruNv+PDhSb8JEyas87WjtDHTBgAAAAAAIIO4aQMAAAAAAJBBdV4epbtYaI7ilGudSs0UxOxq166dZ502aJYe088++8yzTtM1M/v888/X+dxm6TTjWK7xu9/9zjPlUdmmnwtdYd/MbPny5Z7ffvvtOntN5UTPw506dUrahg0b5nnLLbes1vNtt912yeMtttjC8+rVq2vzEpGDTqvu37+/5169eiX9FixY4FnLo/Lxdw866KCkbY899vA8duzYpO2TTz5Z77+Nqml52t1335206ffk4MGDc/ajjHHdDjnkkOTxdddd51mvVy+99NKk38svv+z5+++/L9Cr+5GWaV1zzTVJW6tWrTxrCeuKFSuSfuVaxqg7kB5++OFJm16zjh492nM8v/FbJdv02lN3Ln388ceTfnpNU9Xv1qrGiv6t6v72Pe2005K2KVOmeOZaqjwx0wYAAAAAACCDuGkDAAAAAACQQdy0AQAAAAAAyKA6X9Nm22239XzbbbclbS1atPA8c+bMpO2dd97xPG/ePM+LFi1K+s2fP9/zl19+6TnWaFe31lS39tMcnzNu4VaudP0KXUfBzGz27NmedSvh7777rlrPvXTp0uTxt99+m7OvHhutEaVWP3t0rZS4neJbb73lec2aNXX2mspJs2bNPI8bNy5p22mnnTxrjX+ka5zE7abHjx/vecSIEZ71nI7a0e/M008/3fPcuXOTfvq9mI91FnR7cV0bxcysY8eOnidOnJi0scZD4e29996e4zpUupaCHkO+F6tHz4FxLSel31vPPfdc0pbva0U9prpOjVk6/vTa28zs0Ucf9fzpp596Ltc1bMzS91KP77333pv022STTTx37drV8/HHH5/0y7U+44YU11Mp9XNyrnOemdkOO+zg+Y477vAc1+WLa6wqHS86tuP1qj6O68rpb5unn37a82OPPZb0Yx2bqulvdL12NTMbMmSI53guXLVqleeVK1d6jtcvuh5ZdX+35hszbQAAAAAAADKImzYAAAAAAAAZVCflUTq1TKcvNW/ePOnXpUsXzzrF2szsyCOP9KxTmeLUPp2y9MUXX3jWEgAzs8aNG1f6b8zM1q5d61mn08XyAJ1GNWjQoKTtzTfftHI0Y8YMzzodzczs448/9lybqWW6LaVZOhX4m2++Sdqef/55z6U+/bPYxKmmp5xyiufWrVsnbddff73nQm+PWk60rGbSpEmedaq3WXq+rmqbStWkSZPkcb9+/TwfcMABngcOHJj0i9tsYt0uueQSz23btvUct9qO58f1pVPCY0mGluRUVVKH/Khfv37y+LzzzvMcz7V63K6++urCvrAS1KFDB8+77bZb0qYlRhdccEGl/7tZ/q9HNt10U8/xGmn33Xf3/NVXXyVtWh5VVal5OdFz2Q033OBZy/7N0mO4xx57eB42bFjST7dZ198j8TnqUileD8fznF6DaClbLP3T33B6/aHjPJo2bVryePLkyZ51rC9cuDDpp6Vy8ftYy1P1+JTisVpfep1jZnbZZZd5PvDAA3P2089IfF/12laz/jYxS5dmiedaLYPVMrZ8lx4z0wYAAAAAACCDuGkDAAAAAACQQdy0AQAAAAAAyKB6NamZq1ev3noX2DVo0MBzu3btkjZd70DrdM3M2rdv71lruLt3757005pUrWWMz6eP41oZuv2a/q2mTZsm/bT2La7HMGDAAM/5qEusqKio3oIS65CPY1hdsdZe3+fq1vl17tzZs66XY5ZuMxy3ZNT1dL7++mvPG7hGdHpFRUWPfDxRXR7HfNNxaWY2c+ZMz3HtKV03IG4hv6EU41jU866Z2TPPPONZa/Ljsck1XuL41fUSdE0ws3SNE10jZ9GiRUm/HXfcsdLnK5CiHIvbbLNN8li3oFy+fLnnuB3x4sWL8/o6GjZs6Fm3NzZLj/cuu+yStOV7DBfjWMw3vTYyM5s9e7bnOO712kbXbfjyyy8L8+KqJ7NjMa7J9MADD3jec889k7aRI0d6vvnmmz0XegttPY5PPfVU0qbXw7fcckvSdtVVV3nOx2ssxrGo30dmZoMHD/Z87bXXeq5qTRu9ro3fW8uWLfN89tlnJ20vvPCC53yve7EeMjsWo6233trziBEjkjZdo0k/53FdIf3c6++J+HtRt/KO6z9lcd2ZYhyLcc1EvY649dZbPffv3z/pp78z9Tn0mJml16W6Lq5Z+j2p1zZVvca4/fpLL73k+bTTTsvZrwafl0rHIjNtAAAAAAAAMoibNgAAAAAAABlUJ1t+Ky1VmTNnTtKmj+NUKZ22r1m35DbLvZV3y5Ytk3767z755JOkTbdj0y1VzzzzzKSfTqd79913k7YsTpmra7XZ1tssnZ6mW3fHqcq6hd4VV1yRtOkx5FhkSyyB0mmuugWjWXqMUXu6DbCZWc+ePT3HKeJKp23r9NJ47tbzZJzqrdvf7rPPPp7j50CnwJ5zzjlJWzlv967H54477kja9D3ULSeXLFlS0Ne07bbbeo5bfuu5d8WKFQV9HeVKr49i2UUsiVK6zXMdlCAWJR1ver4yS0ui4nlTy7cLXe6i18CXXnqp5zgWtVTuuuuuS9oKXbZVDOLvgjPOOMOzll3Ea1kdO1qGoVtNm5k1atTI84MPPpi09enTx/OsWbNq8rLLUtyGW0u8W7RokbQ98sgjnrU8parPvB7HeG7kN0Rh6PiIS1z069fPs47FeC2o16ITJkzwPGnSpKSfbscez91bbbWVZ/0+PeaYY5J+rVu39hyXS+nR46dKJj0Pr1mzJum3vp8lZtoAAAAAAABkEDdtAAAAAAAAMqjOy6NUVdOEYpuu2K057nqQa1r4vHnzqv26dOrx+++/X+nfNUt35Ljrrruq/fzlIh5DfV81x12mxo4d61lLZ+Jq4EOHDvUcd0fJ0Gr8COLONrpS/9SpU5M2LbVAzWg54VlnnZW05SqJilNPp0+f7vkvf/mL5yeffDLpp2WvsbRVy09HjRrluVevXkm/Y4891nPc5eHCCy+s9PWWAy1li++ZnhO1TKIQ07n1uJ5wwgmetUzYLN0xhfFbGFoeoztVmqXHKY5nLZFhyn/l9P3r27dv0ta8eXPPsWSmbdu2njfa6Kf/P7Sqkoxc10RR3NXv3HPP9XzUUUd5jteoRx99dM62cqXvZfxe1LJP3Y3voosuSvrpdYqOv2uuuSbpt9lmm3nW3XDMzI444gjPWsbGuPyJXqeMGzcuadPxFq8Xrr76as/VLa3WMcsxKIy4K9fTTz/tOe7Gp+dD/Z0fz8mvvvqqZ/3dV9UxjOdavW+gO3LqGI3/Tj8vZulnVXc9njt3bs7XURvMtAEAAAAAAMggbtoAAAAAAABkEDdtAAAAAAAAMmiDrmmTVbrGyrBhwzzHLadHjhzpma2J101rDLX+b9CgQUm/gw8+2LPWo77yyitJv4kTJ3pm+8ps0/rPuP20rstx2223JW3UFteerr8Qt9fW91XHzvPPP5/0O/nkkz1r3W9N1iNbtGiR55kzZ3ru3bt30m/zzTf3HLcxHj16tOcPP/ww598uBXGNL10vTddIMDO77777POv7XAj6GTrzzDM9xzVt/va3v3lmbbHC2GKLLTzrNqTRggULkseF3gq+FOi6BV26dEna9LpFz1dmZpdddplnPSbxukXXNtFzZby+1LEez4e77LJLpa9X14gwM/vss88MqcaNG3uO6+vp8dDr+wceeCDpp+e122+/3bOeF83MdtxxR89xDYwDDjjA87XXXlvpayh3HTt29Ny9e/ekTb93XnvttaTt448/rvHf4n0vvK5duyaPu3Xr5jmuM6O/C/785z97jsdax2JVa6XqWlbxeninnXbyfOONN3rWNa7M0jEcr22WLl1aac7354qZNgAAAAAAABnETRsAAAAAAIAMojzKfr797fjx4z3rtnJxW+l77rmnsC+shLVs2dLz5ZdfnrTpNDbdcj1uz/jVV18V6NUh33RKeKdOnZI2nbI/Z86cWj1/VdMWy5WWR8XtXhs0aOBZt+QeMmRI0q+6JVFV0Wmuzz77rOdzzjkn6adTW2PpgZbNHXLIIev9mrJs6623Th5vtdVWnuMW2vq+5H0abpjOf+KJJ3pu2LBhztekJVvIHx0fhx9+uOe4HbSe/x599NGkTcciKqflorfeemvS1qFDB896bWiWlkTpNU0cR0q3sp0/f37SpufA9u3bJ236nFpCruWsqNzXX3/tOW5BrI8/+ugjz1WdW/VYxO+tqujnRUt94lby5UbfiwcffNBzPM/pOL333nuTttp8F5bitUTWNGnSJHms565YHqrjQMuUjj/++KRfv379PO+6666eYym5XvPGY63jXv9dLNnSfxeXRLnooos8v/HGG1YozLQBAAAAAADIIG7aAAAAAAAAZBDlUfbzqac6/V4NHTo0ebxmzZpCvaSSpNMeH3roIc+xHECnzOluCFrGYZbN6YxxKnQ5l+ro1ELdISyu6v7kk096jqUW1VXO73MuPXr08BynFuvU0zFjxniOuy7kY4zpc+iU81jeqNNXI53SXurilH0t3125cmXSVsjdYeK57NRTT/Wsx3Tq1KlJv0LvYlWudAwPHjzYczxOq1ev9jxq1KikLYvfmVmj3yXPPPNM0qY7/pxyyilJm+5Astdee3mOJTNaqqpjZ8KECUm/fffd1/Mf/vCHnK/3pptu8qzlVqjc2rVrPc+aNStp02N4/fXXe164cGHST0tzhg8f7llLks3Sa9lIS0z79u3recqUKTn/TTnQ68MddtghZz+9vuzVq1fSpmXYerwjPR9WdW6MZTK5ngNVmzZtWvL4zTff9Bx3B1MDBw70HM+7uZZGqG5Jo1nu4xt/V8yePdvzGWeckbRpSVQhf48w0wYAAAAAACCDuGkDAAAAAACQQdy0AQAAAAAAyCDWtDGzo48+Onmsa6+8/vrrnmPNMWrm/PPP96w137G+UNdtGDdunOdi2K6UtVV+omtxDBgwwHNct+auu+7yTH1w7cVx1Lt3b89xnRRdI0bXVSjEGNN6Yd26MY4VPfZx29PRo0dX2q8UrVixInmsa6fF46hrMOSjplo/Q/r5MTPr2LGjZz0+t99+e9Kv1I/PhtKyZUvPnTt39hzf7xkzZnieN29ewV9XKYvjSNf80m29Ix1Hej0Zn1PXRol0C/C4joOOvyuuuCLnc+DndLxcddVVSdsRRxzhWde6fOGFF3I+h35nvvjii0m/999/3/N+++2XtLVr186zrln01FNPJf2qWhenFOl/7xdffOE5buGs1xVxfRFde0rXzlu8eHHST7cKf+WVVzzHa6lmzZp5jutG6bWU/nYp963bKxPfO92u+7DDDkva+vfv73mPPfbwHK+BdIzNnTvXs64JZma23XbbeY7biyv9fTJ+/Pik7bzzzvMc11msq+seZtoAAAAAAABkEDdtAAAAAAAAMqhsy6N0+8xYHrVq1SrPJ5xwgmdKX2omTmM799xzPevUxjiN8O677/as091q8v7r8zNdf8PYcsstPXfq1Mnz559/nvSL20yjduKUz/33399znKKvY05LKPIxVrQszsxsn3328axlTrrlafzbS5cuTdp0a8hSp1PCzcw++eQTz926dUvaxo4d63nixIme43R+PR/qtP9Ip4WPGDEiadPzuU4h/vDDD3M+H2ovbkO62267edZjEcfsk08+6bkYSopLkV6rVLdMIpZk9OnTx3M8jlpCU27lM/mk15dmaan24MGDPcfvVj0eo0aN8nzhhRfm/Fu77rpr8ljHaY8ePTx37do16VdO331m6ef5zjvv9Dxs2LCkX4MGDTzH3xpaPqrjKp4rDzzwwEr/bvytoduGx+9nvY7R8rizzz476Rf/HdISo4cffjhpmzx5cqX/pqqSUv1dH59Py7sjPb6nn3665wcffDDpl4Xfksy0AQAAAAAAyCBu2gAAAAAAAGQQN20AAAAAAAAyqGzWtIn14Vov16VLl6RNa01Zb6P2Dj300OTxVltt5VnrEidNmpT0u/HGGz1/++23tfrbdVl7yPo5P4pjTI9//fr1Pc+cOTPpx9aI+bH55psnjxs1auQ5rpegdEvM2tLtOC+++OKk7fzzz/esrzGOFa1vvuaaa5K2cqoHj+Ph97//vWddP8Es3UJd12DQ99wsrfVWixYtSh5/9tlnnlu0aFGt1xu3vkR+xDF76qmnVtoWvyOfeOKJwr4wFERc8+Taa6/1HNdx0LUWyvmaY33F906/u3Qr4aOOOirpd91113n+97//7bmqdRfnzJmTPP7ggw88//rXv/Z85plnJv2mT5+e8/WWIv1vvOeeezy3bds26dehQwfPrVu3TtratGnjuarrIP1e1BzfZ13TKK7Fp226ZXz8XtQ1PdesWWNIxfe8qrVrctl+++099+3bN2nTtRbjOL3jjjs8Z/3cykwbAAAAAACADOKmDQAAAAAAQAaVTXnUiSeemDw+5JBDPOv2pWZmZ511lucsTo/KMi2RGT58eNKm2w7/5z//8XzllVcm/XRa4YZ8//W/RXOTJk2Sfvp4/vz5SVttpviVit/85jeedYvMGTNmJP0YY/kR30edChxL17RcTT+/cUqv/jsdv61atUr66Vg/7rjjcv4tfb64je3UqVM9T5gwIWkr58+ITqvv379/0ta9e3fPukVw06ZNk366Peprr73mOde2mmZmf//735PHui2tHjs9vsifWNLWs2fPSvvF7xzdIh7ZpufUSy65JGfb7NmzkzY9VyJ/9Lw2duzYSrNZ7a7rYhnyDTfc4HnMmDGe9beJWVqOs3r16hr/3WL2+eefe47jQ68J4vWBbu98+OGHe44lM/r9qeMtlpyuWLHCcyxD1/O0luDo97FZWrIVS+XK+fpmfTVr1szzyy+/7DluA68lUbNmzUrarrjiCs9ZPxbMtAEAAAAAAMggbtoAAAAAAABkUEmXR22zzTae464bOk3/zjvvTNp0Sh5qRqdy6uruZul7rtP144ruOt1ep6rpvzGruowq145OOn3RLJ3aqOU8Zma//e1vPR900EE5X4dOl7z00kuTtnHjxlm50OmlZmZ77bWXZ11h/9133036ZX06YrFYtWpV8njhwoWed9hhh5z/bsSIEZ51Jwwzs4MPPtizHs84BvSxHmuz3GNx8eLFST/dNYMdiX6i71ncReu5557z/K9//ctzLIfT59BpwlWdN++///6kbbfddlvn60P+xF0tdRq4vuejR49O+tV2x0XUDR1jekz33HPPpJ+WZJxyyilJW212XIzXPvp9rc/HeP5Rvkvb4441L774oufXX3/dcyyr0e/nCy+8sMrnLDX631eTHZf0GlNLkaZNm5b0052C9LdL3DVYl3PQnb7M0msfHTuxbFV3ZmSM1V78naHl3no+jWPj008/9awlc2bFdb3JTBsAAAAAAIAM4qYNAAAAAABABnHTBgAAAAAAIINKbk0brXfTtWri9l9aL6zbfZlRb7g+dD2LuFaNatSokecpU6YkbQsWLPDcuXNnz3HLxLVr13qO6xDp64hbuqsWLVp41jWQzNLPgdahx/8u/W+J25yX05o2cRvo9u3be9b60n/+85919ZLKSqzBHzlypOebbropadN1mE444QTPxx9/fNJPz5v6uY/nyKrOmfq6Vq5c6Tmu07BkyZKcz4F1y8caDHocW7ZsmbTpmhjar5jqwYtJHB/6/ut34YQJE+rsNWH96bVE//79PeuaGmbpOhq1PTfq9XDcQv7777/3zDVv3fvyyy896/phug21mdnAgQM9T5w4MWl78803C/Tqipt+nvX6P66jp1uo6/k1rhml63Pq9X78d3pMb7755qSf/l5Bzei159ChQ5O27bffvtJ/E9f/O+mkkzzr+kLFhpk2AAAAAAAAGcRNGwAAAAAAgAwq+vKoWKqy8847e9577709//DDD0k/3UYvlt2g9nTb4bhtnpY66RThuB1xx44dPVdVYqUlN3F7Ny0V0K0CtSzOzGyLLbao9DWZpdP+X3nlFc9a9mNm1rx5c8+77rprztdb6m699dbksU7HXrRokee4NTUKY8yYMZ6HDBmStLVr185z/fr1PccxEB/n+t9zTUc2M/voo488X3zxxZ6nTp2a8zmwYehU7549eyZtesx1W+mqyk9Re3r9YpaOj+XLl3vWKfnIPv1e1BI4vRYxS8/LTZs2Tdq0zFQ/F3EZAL0OitfA+d7SGjWj16xPPPGE5zjujz32WM+DBg1K2t577z3PVZXf6Lm73L5n9b83Xv9ruVSTJk08d+rUKemnyy3Eax8dRy+99JLnp59+OufrwLrp+6zH44ILLkj66W9ELfm8+uqrk366NXgxn/uYaQMAAAAAAJBB3LQBAAAAAADIIG7aAAAAAAAAZFDRr2mjW9eamV155ZWetb53/vz5Sb977723sC+sTGlN4V577ZW09erVy7PWaPft2zfpp+vC6JaVms3S2vC4vduMGTM8T5o0yXNcT0XXBohrM2ib1sLG9XMaN27sudzWF9A1MOJWlVo3qltO62cEhaM17rqtt5nZX//6V8+6RlODBg2SfrnWlIrrI2ht+F133ZW06TalutViMdcVlyr9zmzVqlXSpsdr6dKlnuP2qMiPOD50TYR58+bl7Ids0+uYtm3beo7n2hYtWngeMGBA0vbYY4951s9FVdc38ZzNGhvZoddEt9xyS9J25JFHej7ssMOStrffftvzfffd5zkeW471j5YtW5Y8vv766yvNbdq0SfrpmI3n21mzZnn+4x//6FnXxETN6Xs+fPhwz3HtL/1sP/XUU57jdaiuw1fMmGkDAAAAAACQQdy0AQAAAAAAyKCiLI/Skow//elPSVufPn0q7Td69OikH1PXCk+3pTQze/zxxyvtN3bs2Lp4OQVTzltYa4maTtU1S7c/jFMVUbdeffXV5PF+++3neeedd/Z80EEHJf122mknzzqe4/M99NBDnuN4oHyjeOh3Ziw51RLR8ePHe+b45o9uc7pw4cKkrUOHDp4XLFjgmdKH4qKlMFrCGrcS1vKAfffdN2nTz8b06dM9xxIA/Vul8DnR81Opnne0hNgs3Tr6mGOOSdp69+7t+f777/ccj7WW3sXy/nISS/OfffZZz1rivf/++yf9Vq9e7fndd99N2t566y3PX331VV5eJ9JSwEMPPdRzLCPV65QhQ4Z4LpVyqIiZNgAAAAAAABnETRsAAAAAAIAMqleTKZP16tXLxPxKXXH/vffeS9oaNmzoWaeeduvWLek3d+7cAr26wqioqKi37l7rlpVjWKamV1RU9MjHE2XxOMbp3aoUpmb/P8ZiSSjpsVhbWpKx5557Jm1bb72153/84x+eN+SOeaU8Frt06ZI8vuiiizzr7kFTpkxJ+hVh2UhZjUX9njz99NM9X3LJJUk/LXvSnW3MzD744APPWuofS1/0e7fQ38GlPBY3pB133NGzlvOYmb3zzjue9bO0ZMmS2v65shqLpaoYx+Jmm22WPJ42bZrnXXbZJee/e/jhhz0PGjTIc9wtrwhVOhaZaQMAAAAAAJBB3LQBAAAAAADIIG7aAAAAAAAAZFDRrGmjW/2NGzfO83HHHZfz33z00Ueeu3btmrQV29ZsxVijiJ+hXrgEMBZLAmOxhnQtjqysUcVYLAllOxZLaR04xmJh6Gdk8uTJSduCBQs863bH66Fsx2IpKcaxqGvmmZnNnDnTc7NmzTyvWrUq6bf77rt7/lgcRSAAAAF3SURBVPDDDwv06jYI1rQBAAAAAAAoFty0AQAAAAAAyKCN190lG3QrQ92u+5tvvkn6rVixwvMFF1zg+bvvvivgqwMAoHQVW7kGkHWMKayLfkYGDhyYtPG7BqVi2bJlyePx48d7HjBggOeTTz456afLoJQDZtoAAAAAAABkEDdtAAAAAAAAMoibNgAAAAAAABlUNFt+q402+uleU9wyUde+KaV64WLcwg0/w3aKJYCxWBIYiyWAsVgSGIslgLFYEhiLJYCxWBLY8hsAAAAAAKBYcNMGAAAAAAAgg2q65fdyM5tfiBdSE1oCVSba5fG5MnEMyxTHsfhxDEsDx7H4cQxLA8ex+HEMSwPHsfhxDEtDpcexRmvaAAAAAAAAoG5QHgUAAAAAAJBB3LQBAAAAAADIIG7aAAAAAAAAZBA3bQAAAAAAADKImzYAAAAAAAAZxE0bAAAAAACADOKmDQAAAAAAQAZx0wYAAAAAACCDuGkDAAAAAACQQf8HAb3G5WzhOvsAAAAASUVORK5CYII=\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: 4.12 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(8*1,8)) # 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.0018, PSNR (dB) = 27.48\n",
      "MSE = 0.0040, PSNR (dB) = 23.95\n",
      "MSE = 0.0007, PSNR (dB) = 31.78\n",
      "MSE = 0.0031, PSNR (dB) = 25.14\n",
      "MSE = 0.0037, PSNR (dB) = 24.35\n",
      "MSE = 0.0005, PSNR (dB) = 33.02\n",
      "MSE = 0.0046, PSNR (dB) = 23.34\n",
      "MSE = 0.0042, PSNR (dB) = 23.81\n",
      "MSE = 0.0061, PSNR (dB) = 22.12\n",
      "MSE = 0.0044, PSNR (dB) = 23.60\n",
      "Wall time: 3.01 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.9815\n",
      "0.9731\n",
      "0.9896\n",
      "0.9738\n",
      "0.9498\n",
      "0.9927\n",
      "0.9296\n",
      "0.9573\n",
      "0.9324\n",
      "0.9617\n",
      "Wall time: 369 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.9966\n",
      "0.9951\n",
      "0.9979\n",
      "0.9948\n",
      "0.9901\n",
      "0.9990\n",
      "0.9908\n",
      "0.9928\n",
      "0.9834\n",
      "0.9880\n",
      "Wall time: 943 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": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Deep Autoencoders in Keras - 784-256-128-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
}