Deep Autoencoders in Keras - 784-256-128-64-128-256-784-tf2-v2-reuse.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# RPi is very memory bound, import only what we need\n",
    "from tensorflow.keras.models import Model, save_model, load_model\n",
    "from tensorflow.keras.layers import Input, Dense\n",
    "from tensorflow.keras.datasets import mnist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "MODEL_VERSION = 'v2'\n",
    "INPUT_SIZE = (784,)\n",
    "BOTTLENECK_DIM = 64\n",
    "BOTTLENECK_SIZE = (BOTTLENECK_DIM,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_train shape = (60000, 784)\n",
      "x_test shape = (10000, 784)\n",
      "Wall time: 584 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import numpy as np\n",
    "\n",
    "# Load mnist dataset\n",
    "(x_train, y_train),(x_test, y_test) = mnist.load_data()\n",
    "\n",
    "# Preprocess\n",
    "x_train = x_train.astype('float32')/255 # normalise to 0 and 1\n",
    "x_test  = x_test.astype('float32')/255\n",
    "x_train = x_train.reshape(len(x_train), (np.prod(x_train.shape[1:]))) # 28x28 flatten into 1x784\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}') # 60000 images 784 vectors = 60000 row 784 col\n",
    "print(f'x_test shape = {x_test.shape}') # 10000 images 784 vectors = 10000 row 784 col"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:No training configuration found in the save file, so the model was *not* compiled. Compile it manually.\n",
      "WARNING:tensorflow:No training configuration found in the save file, so the model was *not* compiled. Compile it manually.\n",
      "Wall time: 266 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# load model\n",
    "# Error message see\n",
    "# https://stackoverflow.com/questions/53295570/userwarning-no-training-configuration-found-in-save-file-the-model-was-not-c\n",
    "#autoencoder = load_model(f'autoencoder{MODEL_VERSION}/model/autoencoder{MODEL_VERSION}-784-256-128-64-128-256-784-best.h5')\n",
    "encoder = load_model(f'autoencoder{MODEL_VERSION}/model/encoder{MODEL_VERSION}-784-256-128-64-best.h5')\n",
    "decoder = load_model(f'autoencoder{MODEL_VERSION}/model/decoder{MODEL_VERSION}-64-128-256-784-best.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"encoder\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_2 (InputLayer)         [(None, 784)]             0         \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 256)               200960    \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 128)               32896     \n",
      "_________________________________________________________________\n",
      "dense_3 (Dense)              (None, 64)                8256      \n",
      "=================================================================\n",
      "Total params: 242,112\n",
      "Trainable params: 242,112\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n",
      "Model: \"decoder\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_3 (InputLayer)         [(None, 64)]              0         \n",
      "_________________________________________________________________\n",
      "dense_4 (Dense)              (None, 128)               8320      \n",
      "_________________________________________________________________\n",
      "dense_5 (Dense)              (None, 256)               33024     \n",
      "_________________________________________________________________\n",
      "dense_6 (Dense)              (None, 784)               201488    \n",
      "=================================================================\n",
      "Total params: 242,832\n",
      "Trainable params: 242,832\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "#autoencoder.summary()\n",
    "encoder.summary()\n",
    "decoder.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 1.41 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "encoded_imgs = encoder.predict(x_test)\n",
    "decoded_imgs = decoder.predict(encoded_imgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(784,)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAOBUlEQVR4nO3dfahcdX7H8c/HmChkI6jBGI2p20WxpWC2BC3Eh5RlJSpi9o+Y9bl04W50lY2INmz/UCiV0FbrH4Lkho2bFpvNgg8ry2IiQZuKuOSBVOPqrmmwmuSaaEPcLBptkm//uOdur/HOb+6dpzM33/cLhpk53zkz3wz3k3NmfnPOzxEhACe/U+puAEBvEHYgCcIOJEHYgSQIO5DEqb18Mdt89Q90WUR4rOVtbdltL7L9G9u7bK9o57kAdJdbHWe3PUXSbyV9W9IeSVsk3RwRvy6sw5Yd6LJubNkvk7QrInZHxBeSfirpxjaeD0AXtRP28yV9MOr+nmrZl9gesL3V9tY2XgtAm9r5gm6sXYWv7KZHxKCkQYndeKBO7WzZ90i6YNT9OZL2tdcOgG5pJ+xbJF1k++u2p0n6rqQXOtMWgE5reTc+Io7avkfSBklTJK2JiLc61hmAjmp56K2lF+MzO9B1XflRDYDJg7ADSRB2IAnCDiRB2IEkCDuQRE+PZ0fvPfTQQ8X6HXfcUawvXbq0WN+6lUMeJgu27EAShB1IgrADSRB2IAnCDiRB2IEkGHo7CSxcuLBhbWBgoLjup59+WqzPnz+/WGfobfJgyw4kQdiBJAg7kARhB5Ig7EAShB1IgrADSXB22UlgxowZxfru3bsb1tauXVtcd8WK8uS7zf4+jh07Vqyj9zi7LJAcYQeSIOxAEoQdSIKwA0kQdiAJwg4kwfHsk8Bdd91VrB85cqRh7dFHHy2ue/To0ZZ6wuTTVthtvyfpsKRjko5GRPlMBwBq04kt+19GxMcdeB4AXcRndiCJdsMekjba3mZ7zJOd2R6wvdU2JysDatTubvyCiNhn+xxJL9l+JyI2j35ARAxKGpQ4EAaoU1tb9ojYV10fkPScpMs60RSAzms57Lan254xclvSNZJ2dqoxAJ3Vzm78LEnP2R55nn+LiBc70hW+5MEHHyzWV61a1bA2NDTU6XYwSbUc9ojYLenSDvYCoIsYegOSIOxAEoQdSIKwA0kQdiAJDnHtA81OFX3aaacV6++8804n28FJii07kARhB5Ig7EAShB1IgrADSRB2IAnCDiTBOHsfWLRoUVvrv/giRxajObbsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AE4+x9YNmyZcX6559/Xqx/9NFHnWwHJym27EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBOPsPVBNa93Q2WefXaxv2rSpk+30jYULFxbrS5cubev5Dx061LC2efPm4rrNzhEQES31VKemW3bba2wfsL1z1LKzbL9k+93q+szutgmgXePZjf+JpBNPpbJC0qaIuEjSpuo+gD7WNOwRsVnSwRMW3yhpbXV7raTFHe4LQIe1+pl9VkQMSVJEDNk+p9EDbQ9IGmjxdQB0SNe/oIuIQUmDkmR78n2rAZwkWh162297tiRV1wc61xKAbmg17C9IurO6faekn3emHQDd4mbjhbbXSVooaaak/ZIekvS8pJ9JmivpfUlLIuLEL/HGeq6Uu/HnnXdesb5nz55i/dZbby3W161bN+GeOmXatGnF+sqVKxvWli9fXlz3/fffL9YPHz7c8vpXXHFFcd0lS5YU6xs3bizW6xQRY/6wo+ln9oi4uUHpW211BKCn+LkskARhB5Ig7EAShB1IgrADSXCI6yRQ56miTzmlvD1YvXp1sX777bc3rN19993FdZ966qlivdkptksWLy4fzrFq1apifd68ecX6J598MuGeuo0tO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kwTh7D8ydO7et9bds2dKhTibuiSeeKNavueaaluvNTpHdzdM1b9iwoVg//fTTi/Xp06cX64yzA6gNYQeSIOxAEoQdSIKwA0kQdiAJwg4kwTh7D8yaNavuFho699xzi/UbbrihWL/llluK9ZdffnnCPfXCZ599Vqzv2rWrWL/yyiuL9fXr10+4p25jyw4kQdiBJAg7kARhB5Ig7EAShB1IgrADSTDO3gNffPFFW+vPmTOnWG/n2OnbbrutWG82Dv/aa6+1/NqT2YwZM+puYcKabtltr7F9wPbOUcsetr3X9o7qcl132wTQrvHsxv9E0qIxlv9zRMyrLr/sbFsAOq1p2CNis6SDPegFQBe18wXdPbbfqHbzz2z0INsDtrfa3trGawFoU6thf1LSNyTNkzQk6dFGD4yIwYiYHxHzW3wtAB3QUtgjYn9EHIuI45JWS7qss20B6LSWwm579qi735G0s9FjAfSHpuPsttdJWihppu09kh6StND2PEkh6T1J3+9ij5Peq6++Wqx/+OGHxfqyZcuK9XvvvXfCPY14/fXXi/VTTy3/iVx99dXF+saNGyfcUy80+3edccYZxfqhQ4c62U5PNA17RNw8xuIfd6EXAF3Ez2WBJAg7kARhB5Ig7EAShB1IgkNce+Dw4cPF+t69e4v1JUuWFOv33Xdfw9rRo0eL6x48WD7s4fjx48X6lClTivV+1Wy4stmhvc2mm+5HbNmBJAg7kARhB5Ig7EAShB1IgrADSRB2IAlHRO9ezO7di00iS5cuLdaffvrpYv3JJ59sWGvn8FdJGhwcLNavv/76Yn3NmjUNa0eOHGmppxHNDh2eO3duw9rq1auL61577bXFer9ORS1JEeGxlrNlB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkGGefBNavX1+sL168uGHt8ccfL6772GOPFevNpoNetGisOT//38yZMxvW7DGHg/9g2rRpxfrFF19crF966aUNa/fff39x3W3bthXr/YxxdiA5wg4kQdiBJAg7kARhB5Ig7EAShB1IgnH2SWDq1KnF+iOPPNKwtnz58uK6zc5Z//zzzxfrH3zwQbFeUvp9gCQtWLCgWG927vYHHnigYW3Hjh3FdSezlsfZbV9g+2Xbb9t+y/YPq+Vn2X7J9rvV9ZmdbhpA54xnN/6opPsj4k8k/YWkH9j+U0krJG2KiIskbaruA+hTTcMeEUMRsb26fVjS25LOl3SjpLXVw9ZKKu+TAajVhOZ6s32hpG9K+pWkWRExJA3/h2D7nAbrDEgaaK9NAO0ad9htf03SM5KWR8Tvmh3EMCIiBiUNVs/BF3RATcY19GZ7qoaD/nREPFst3m97dlWfLelAd1oE0AlNh948vAlfK+lgRCwftfwfJf1PRKy0vULSWRHxYJPnYsveY5dffnmxftNNNxXrV111VbF+ySWXFOuvvPJKw9r27duL627evLlYb3Y652bTTZ+sGg29jWc3foGk2yW9aXtkcPJHklZK+pnt70l6X1J5EnEAtWoa9oh4VVKjD+jf6mw7ALqFn8sCSRB2IAnCDiRB2IEkCDuQBIe4AicZTiUNJEfYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJNA277Qtsv2z7bdtv2f5htfxh23tt76gu13W/XQCtajpJhO3ZkmZHxHbbMyRtk7RY0k2Sfh8R/zTuF2OSCKDrGk0SMZ752YckDVW3D9t+W9L5nW0PQLdN6DO77QslfVPSr6pF99h+w/Ya22c2WGfA9lbbW9vqFEBbxj3Xm+2vSfp3SX8fEc/aniXpY0kh6e80vKv/102eg914oMsa7caPK+y2p0r6haQNEfHYGPULJf0iIv6syfMQdqDLWp7Y0bYl/VjS26ODXn1xN+I7kna22ySA7hnPt/FXSPoPSW9KOl4t/pGkmyXN0/Bu/HuSvl99mVd6LrbsQJe1tRvfKYQd6D7mZweSI+xAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiTR9ISTHfaxpP8edX9mtawf9Wtv/dqXRG+t6mRvf9So0NPj2b/y4vbWiJhfWwMF/dpbv/Yl0VuretUbu/FAEoQdSKLusA/W/Pol/dpbv/Yl0VuretJbrZ/ZAfRO3Vt2AD1C2IEkagm77UW2f2N7l+0VdfTQiO33bL9ZTUNd6/x01Rx6B2zvHLXsLNsv2X63uh5zjr2aeuuLabwL04zX+t7VPf15zz+z254i6beSvi1pj6Qtkm6OiF/3tJEGbL8naX5E1P4DDNtXSfq9pH8ZmVrL9j9IOhgRK6v/KM+MiL/pk94e1gSn8e5Sb42mGf8r1fjedXL681bUsWW/TNKuiNgdEV9I+qmkG2voo+9FxGZJB09YfKOktdXttRr+Y+m5Br31hYgYiojt1e3DkkamGa/1vSv01RN1hP18SR+Mur9H/TXfe0jaaHub7YG6mxnDrJFptqrrc2ru50RNp/HupROmGe+b966V6c/bVUfYx5qapp/G/xZExJ9LulbSD6rdVYzPk5K+oeE5AIckPVpnM9U0489IWh4Rv6uzl9HG6Ksn71sdYd8j6YJR9+dI2ldDH2OKiH3V9QFJz2n4Y0c/2T8yg251faDmfv4gIvZHxLGIOC5ptWp876ppxp+R9HREPFstrv29G6uvXr1vdYR9i6SLbH/d9jRJ35X0Qg19fIXt6dUXJ7I9XdI16r+pqF+QdGd1+05JP6+xly/pl2m8G00zrprfu9qnP4+Inl8kXafhb+T/S9Lf1tFDg77+WNJ/Vpe36u5N0joN79b9r4b3iL4n6WxJmyS9W12f1Ue9/auGp/Z+Q8PBml1Tb1do+KPhG5J2VJfr6n7vCn315H3j57JAEvyCDkiCsANJEHYgCcIOJEHYgSQIO5AEYQeS+D/Y7GpN55lJcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(64,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAALzElEQVR4nO3dX4hc9RnG8ecxf9yYZI20toorjUIISLFRlqAElGpbtEpSoYIBxWrRK0VJQUxzIb1RUBF7UYQkakWt0q6KIVitEMUK1mr+tI1uLDZY3Bobg0pMmjRG317sRNZm456ZOec3k5fvBxZ3doZ5nzF5cs6cPXN+jggByOOYXgcAUC9KDSRDqYFkKDWQDKUGkpnexJMODAzEnDlzmnjqwxx77LFF5kjSvn37is0q7aOPPio2a/r0Rv7aTWrBggXFZknSnj17isz58MMPtWfPHk92XyP/d+fMmaOlS5c28dSHmT9/fpE5kjQ6OlpsliQdPHiw2KyRkZFis0444YRisx5++OFisyTplVdeKTLnzjvvPOJ97H4DyVBqIBlKDSRDqYFkKDWQDKUGkqHUQDKUGkiGUgPJVCq17Ytsv2X7bdu3Nh0KQOemLLXtaZJ+JeliSWdIWm77jKaDAehMlS31YklvR8T2iDgg6XFJy5qNBaBTVUp9iqR3J9wea/3sS2xfb/t126/v37+/rnwA2lSl1JN9vOuwqxVGxOqIGI6I4YGBge6TAehIlVKPSTp1wu0hSe81EwdAt6qU+jVJC2yfZnumpCskrWs2FoBOTXmRhIg4aPsGSc9JmibpgYh4o/FkADpS6conEfGMpGcazgKgBpxRBiRDqYFkKDWQDKUGkqHUQDKUGkiGUgPJNLJCx7x584qt0HHZZZcVmSNJt99+e7FZklTyHPqSH8JZv359sVlr1qwpNkuSli9fXmTOrFmzjngfW2ogGUoNJEOpgWQoNZAMpQaSodRAMpQaSIZSA8lQaiAZSg0kU2WFjgds77S9tUQgAN2psqX+taSLGs4BoCZTljoiXpL0YYEsAGpQ23vqicvu7N69u66nBdCm2ko9cdmdwcHBup4WQJs4+g0kQ6mBZKr8SusxSa9IWmh7zPZPm48FoFNV1tIqc30WALVg9xtIhlIDyVBqIBlKDSRDqYFkKDWQDKUGkmlk2Z19+/Zp69Z8H79euXJl0Xm2i80655xzis0q6brrris6b3h4uOi8ybClBpKh1EAylBpIhlIDyVBqIBlKDSRDqYFkKDWQDKUGkqHUQDJVrlF2qu0XbI/afsP2TSWCAehMlXO/D0r6WURssj1X0kbbz0fEmw1nA9CBKsvu7IiITa3vP5E0KumUpoMB6Exb76ltz5d0lqRXJ7nvi2V39u7dW086AG2rXGrbcyQ9IenmiDhssayJy+7Mnj27zowA2lCp1LZnaLzQj0bEk81GAtCNKke/Lel+SaMRcU/zkQB0o8qWeomkqyRdYHtL6+uHDecC0KEqy+68LKncdXUAdIUzyoBkKDWQDKUGkqHUQDKUGkiGUgPJUGogGUoNJOOIqP9J7fqftA+sWrWq6LyBgYFisxYtWlRs1v79+4vNuvzyy4vNkqSrr766yJz169dr165dk54UxpYaSIZSA8lQaiAZSg0kQ6mBZCg1kAylBpKh1EAylBpIpsqFBwds/9n2X1rL7vyiRDAAnamy7M5/JV0QEXtalwp+2fbvI+JPDWcD0IEqFx4MSXtaN2e0vlKe2w1kUPVi/tNsb5G0U9LzEfGVy+7UHRJAdZVKHRGfRcQiSUOSFtv+9iSP+WLZnbpDAqiuraPfEfGxpBclXdRIGgBdq3L0+0Tb81rfz5L0PUnbmg4GoDNVjn6fLOkh29M0/o/AbyNifbOxAHSqytHvv2p8TWoARwHOKAOSodRAMpQaSIZSA8lQaiAZSg0kQ6mBZCg1kEyVM8r62rJly4rNOv7444vNkqTdu3cXm3X33XcXmzU0NFRs1syZM4vNkqTh4TKfZ9qwYcMR72NLDSRDqYFkKDWQDKUGkqHUQDKUGkiGUgPJUGogGUoNJEOpgWQql7p1Qf/NtrnoINDH2tlS3yRptKkgAOpRddmdIUmXSFrbbBwA3aq6pb5X0i2SPj/SA1hLC+gPVVbouFTSzojY+FWPYy0toD9U2VIvkbTU9juSHpd0ge1HGk0FoGNTljoiVkbEUETMl3SFpA0RcWXjyQB0hN9TA8m0dTmjiHhR40vZAuhTbKmBZCg1kAylBpKh1EAylBpIhlIDyVBqIJlGlt2ZO3euFi9e3MRTH+auu+4qMkeSVqxYUWyWJJ1//vnFZl177bXFZu3cubPYrMHBwWKzJGnBggVF5gwMDBzxPrbUQDKUGkiGUgPJUGogGUoNJEOpgWQoNZAMpQaSodRAMpQaSKbSaaKtK4l+IukzSQe5DDDQv9o59/u7EbGrsSQAasHuN5BM1VKHpD/Y3mj7+skeMHHZnU8//bS+hADaUnX3e0lEvGf7G5Ket70tIl6a+ICIWC1ptSQNDg5GzTkBVFRpSx0R77X+u1PSU5LKfFgaQNuqLJA32/bcQ99L+oGkrU0HA9CZKrvf35T0lO1Dj/9NRDzbaCoAHZuy1BGxXdJ3CmQBUAN+pQUkQ6mBZCg1kAylBpKh1EAylBpIhlIDyTSy7I4kHXNMmX8vbrvttiJzJGnhwoXFZknSmjVris0688wzi81qnchUxObNm4vNkqSTTjqpyJy9e/ce8T621EAylBpIhlIDyVBqIBlKDSRDqYFkKDWQDKUGkqHUQDKUGkimUqltz7M9Ynub7VHb5zYdDEBnqp77/UtJz0bEj23PlHRcg5kAdGHKUtselHSepJ9IUkQckHSg2VgAOlVl9/t0SR9IetD2ZttrW9f//hKW3QH6Q5VST5d0tqT7IuIsSXsl3fr/D4qI1RExHBHDM2bMqDkmgKqqlHpM0lhEvNq6PaLxkgPoQ1OWOiLel/Su7UNXCLhQ0puNpgLQsapHv2+U9GjryPd2Sdc0FwlANyqVOiK2SBpuOAuAGnBGGZAMpQaSodRAMpQaSIZSA8lQaiAZSg0kQ6mBZBpbSyujp59+utcRGjMyMlJs1h133FFs1o4dO4rNkqRVq1YVmbNu3boj3seWGkiGUgPJUGogGUoNJEOpgWQoNZAMpQaSodRAMpQaSGbKUtteaHvLhK/dtm8uEQ5A+6Y8TTQi3pK0SJJsT5P0L0lPNZwLQIfa3f2+UNI/IuKfTYQB0L12S32FpMcmu4Nld4D+ULnUrWt+L5X0u8nuZ9kdoD+0s6W+WNKmiPh3U2EAdK+dUi/XEXa9AfSPSqW2fZyk70t6stk4ALpVddmd/0j6WsNZANSAM8qAZCg1kAylBpKh1EAylBpIhlIDyVBqIBlKDSTjiKj/Se0PJLX78cyvS9pVe5j+kPW18bp651sRceJkdzRS6k7Yfj0ihnudowlZXxuvqz+x+w0kQ6mBZPqp1Kt7HaBBWV8br6sP9c17agD16KctNYAaUGogmb4ote2LbL9l+23bt/Y6Tx1sn2r7Bdujtt+wfVOvM9XJ9jTbm22v73WWOtmeZ3vE9rbWn925vc7Urp6/p24tEPB3jV8uaUzSa5KWR8SbPQ3WJdsnSzo5IjbZnitpo6QfHe2v6xDbKyQNSxqMiEt7nacuth+S9MeIWNu6gu5xEfFxr3O1ox+21IslvR0R2yPigKTHJS3rcaauRcSOiNjU+v4TSaOSTultqnrYHpJ0iaS1vc5SJ9uDks6TdL8kRcSBo63QUn+U+hRJ7064PaYkf/kPsT1f0lmSXu1tktrcK+kWSZ/3OkjNTpf0gaQHW28t1tqe3etQ7eqHUnuSn6X5PZvtOZKekHRzROzudZ5u2b5U0s6I2NjrLA2YLulsSfdFxFmS9ko66o7x9EOpxySdOuH2kKT3epSlVrZnaLzQj0ZElssrL5G01PY7Gn+rdIHtR3obqTZjksYi4tAe1YjGS35U6YdSvyZpge3TWgcmrpC0rseZumbbGn9vNhoR9/Q6T10iYmVEDEXEfI3/WW2IiCt7HKsWEfG+pHdtL2z96EJJR92BzUrX/W5SRBy0fYOk5yRNk/RARLzR41h1WCLpKkl/s72l9bOfR8QzPcyEqd0o6dHWBma7pGt6nKdtPf+VFoB69cPuN4AaUWogGUoNJEOpgWQoNZAMpQaSodRAMv8DHf/I9wrUCscAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(784,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAO+klEQVR4nO3df4xV9ZnH8c8DFDEUlFEEAkQQIVltolUgayQKaaiuiSKYmvKH0dQEEiEpatwl+AdE3MTsLup/TSAl4KYr1qgplipVUnT9BxkUEDpSkbDDlAFkUaFK5Nezf8xhM8U53zPcX+cOz/uVTO6957nn3ofLfOace7/3nK+5uwBc+vqV3QCAxiDsQBCEHQiCsANBEHYgiAGNfDIz46N/oM7c3XpaXtWW3czuNrM9ZrbXzBZX81gA6ssqHWc3s/6S/iJppqQOSVslzXX3PyfWYcsO1Fk9tuxTJe11933ufkrSOkmzqng8AHVUTdhHSzrQ7XZHtuzvmNk8M2s1s9YqngtAlar5gK6nXYXv7aa7+0pJKyV244EyVbNl75A0ttvtMZIOVtcOgHqpJuxbJU00s/FmNlDSzyWtr01bAGqt4t14dz9jZgslbZTUX9Jqd99ds84A1FTFQ28VPRnv2YG6q8uXagD0HYQdCIKwA0EQdiAIwg4EQdiBIBp6PDvq4/LLL8+tPfnkk8l129vbk/UPP/wwWd+7d2+yfubMmWQdjcOWHQiCsANBEHYgCMIOBEHYgSAIOxAEQ299QEtLS7I+Z86c3Nr8+fOT665bty5ZX78+fYoChtb6DrbsQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAE4+x9wNixY5P1Bx54ILe2YcOG5LrLly9P1o8fP56so+9gyw4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQTDO3gT69Uv/zX3kkUeSdbMeJ+2UJL3yyivJdU+cOJGs49JRVdjNbL+kE5LOSjrj7pNr0RSA2qvFln2Gux+tweMAqCPeswNBVBt2l/RHM9tmZvN6uoOZzTOzVjNrrfK5AFSh2t342939oJldI+kdM/vU3d/vfgd3XylppSSZmVf5fAAqVNWW3d0PZpdHJL0haWotmgJQexWH3cwGm9mQ89cl/VTSrlo1BqC2qtmNHyHpjWyMd4Ck/3L3t2vSVTDDhw9P1m+88cZkfefOnbm1LVu2JNd1551VFBWH3d33Sbqphr0AqCOG3oAgCDsQBGEHgiDsQBCEHQiCQ1ybwIwZM5L1oUOHJuubN2/OrZ08ebKSlhqi6NDec+fONaiTGNiyA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQjLM3QOpUz5I0ZcqUZH306NHJ+qeffppb68uHsBa9bn3531YGtuxAEIQdCIKwA0EQdiAIwg4EQdiBIAg7EATj7E1gwoQJyfrp06eT9QMHDtSynYYpGie/7LLLkvUBA9K/vt98803Fz30pYssOBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0Ewzt4AgwYNStbHjx+frG/dujVZP3Xq1EX31Cj9+/fPrd17773JdZ966qlkfeTIkcn6jh07cmsrVqxIrls01fWZM2eS9WZUuGU3s9VmdsTMdnVb1mJm75jZZ9nlsPq2CaBavdmNXyPp7guWLZa0yd0nStqU3QbQxArD7u7vSzp2weJZktZm19dKur/GfQGosUrfs49w905JcvdOM7sm745mNk/SvAqfB0CN1P0DOndfKWmlJJlZvKMPgCZR6dDbYTMbJUnZ5ZHatQSgHioN+3pJD2fXH5b0u9q0A6BeCnfjzexlSdMlXW1mHZKWSnpO0m/N7FFJ7ZJ+Vs8m+7prr702WR8yZEiyvn379mS9zGOzx40bl6yvWrUqtzZ16tTkut99912yfvTo0WT9yiuvzK3Nnj07uW7ROQRaW1uT9WacW74w7O4+N6f0kxr3AqCO+LosEARhB4Ig7EAQhB0IgrADQXCIawNMmjQpWS86ZfKXX36ZrKemNq52WG7s2LHJ+saNG5P1UaNG5dbWrl2bW5OkpUuXJuupU0VL0pgxY3Jry5YtS6570003JettbW3J+okTJ5L1MrBlB4Ig7EAQhB0IgrADQRB2IAjCDgRB2IEgGGdvgKFDhybrRaclPnToULJezVh6v37pv/dz5sxJ1ltaWpL1F154Ibf2zDPPJNc9e/Zssl5k3759ubVXX301ue6MGTOS9dtuuy1Zf/fdd5P1Mg6BZcsOBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0Ewzt4A7e3tyXrROHk9j2cfPnx4sv7QQw8l62+++Wayvnz58txatePoRVL/9t27dyfXve+++5L1O++8M1nfvHlzsl7GNNts2YEgCDsQBGEHgiDsQBCEHQiCsANBEHYgCMbZG+D48ePJetF4c9Gx0x988EFurei46YULFybrqXOvS9Jjjz2WrBcdq1+Wr7/+Olnv7OxM1ou+v9CMUzYXbtnNbLWZHTGzXd2WLTOzv5rZ9uznnvq2CaBavdmNXyPp7h6Wv+DuN2c/f6htWwBqrTDs7v6+pGMN6AVAHVXzAd1CM9uZ7eYPy7uTmc0zs1Yza63iuQBUqdKw/0rSBEk3S+qUtCLvju6+0t0nu/vkCp8LQA1UFHZ3P+zuZ939nKRVkqbWti0AtVZR2M2s+zy8syXtyrsvgOZQOM5uZi9Lmi7pajPrkLRU0nQzu1mSS9ovaX4de+zzOjo6kvUhQ4Yk67NmzUrWX3zxxdza6dOnk+teccUVyXrRdwBuueWWZH3btm0VP3Y9FZ0v//PPP6/q8cv8t+UpDLu7z+1h8a/r0AuAOuLrskAQhB0IgrADQRB2IAjCDgTBIa4N8NVXXyXrb7/9drI+dWr6O0tXXXVVbq3oUM3nn38+WZ8yZUqyvmjRomT9448/zq21tqa/QV3t8NXAgQNza9OnT0+ue+uttybra9asSdarmUa7XtiyA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQjLM3QNHplF966aVk/a677krWn3jiidza008/nVz30KFDyfqzzz6brD/++OPJempq46JDf0+ePJmsDx48OFm/4447cmtFU1EXjaPv2LEjWW9GbNmBIAg7EARhB4Ig7EAQhB0IgrADQRB2IAhr5HG3ZtZ8B/k2gQED0l932LNnT7J+3XXX5daWLFmSXHf16tXJupkl69dff32ynvqOwLfffptcd+LEicn6DTfckKwfPnw4t5aa5lpKn55bKj5Fd5ncvcf/NLbsQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAE4+x9QEtLS7KeOv/6iBEjkuu+9dZbyfqGDRuS9WnTpiXrM2fOzK0VTRf9xRdfJOttbW3J+oIFC3Jr7e3tyXX7sorH2c1srJn9yczazGy3mf0yW95iZu+Y2WfZ5bBaNw2gdnqzG39G0pPu/g+S/lHSAjO7QdJiSZvcfaKkTdltAE2qMOzu3unuH2XXT0hqkzRa0ixJa7O7rZV0f72aBFC9izoHnZmNk/RjSVskjXD3TqnrD4KZXZOzzjxJ86prE0C1eh12M/uhpNckLXL340UHSJzn7islrcwegw/ogJL0aujNzH6grqD/xt1fzxYfNrNRWX2UpCP1aRFALRQOvVnXJnytpGPuvqjb8n+X9L/u/pyZLZbU4u7/XPBYbNnrYNCgQbm1Bx98MLnuyJEjk/X+/fsn65MmTUrWU4eCFu0drlu3Lll/7733kvWiU3hfqvKG3nqzG3+7pIckfWJm27NlSyQ9J+m3ZvaopHZJP6tFowDqozDs7v6BpLw/wT+pbTsA6oWvywJBEHYgCMIOBEHYgSAIOxAEh7giqbfflKxEI3/3IuFU0kBwhB0IgrADQRB2IAjCDgRB2IEgCDsQxEWdlgrxMBZ+6WDLDgRB2IEgCDsQBGEHgiDsQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0EUht3MxprZn8yszcx2m9kvs+XLzOyvZrY9+7mn/u0CqFThJBFmNkrSKHf/yMyGSNom6X5JD0r6m7v/R6+fjEkigLrLmySiN/Ozd0rqzK6fMLM2SaNr2x6Aeruo9+xmNk7SjyVtyRYtNLOdZrbazIblrDPPzFrNrLWqTgFUpddzvZnZDyW9J+lf3f11Mxsh6agkl7RcXbv6vyh4DHbjgTrL243vVdjN7AeSfi9po7s/30N9nKTfu/uPCh6HsAN1VvHEjtY1jeevJbV1D3r2wd15syXtqrZJAPXTm0/jp0n6b0mfSDqXLV4iaa6km9W1G79f0vzsw7zUY7FlB+qsqt34WiHsQP0xPzsQHGEHgiDsQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0EQdiCIwhNO1thRSf/T7fbV2bJm1Ky9NWtfEr1Vqpa9XZtXaOjx7N97crNWd59cWgMJzdpbs/Yl0VulGtUbu/FAEIQdCKLssK8s+flTmrW3Zu1LordKNaS3Ut+zA2icsrfsABqEsANBlBJ2M7vbzPaY2V4zW1xGD3nMbL+ZfZJNQ13q/HTZHHpHzGxXt2UtZvaOmX2WXfY4x15JvTXFNN6JacZLfe3Knv684e/Zzay/pL9ImimpQ9JWSXPd/c8NbSSHme2XNNndS/8ChpndIelvkl46P7WWmf2bpGPu/lz2h3KYu/9Lk/S2TBc5jXedesubZvwRlfja1XL680qUsWWfKmmvu+9z91OS1kmaVUIfTc/d35d07ILFsyStza6vVdcvS8Pl9NYU3L3T3T/Krp+QdH6a8VJfu0RfDVFG2EdLOtDtdoeaa753l/RHM9tmZvPKbqYHI85Ps5VdXlNyPxcqnMa7kS6YZrxpXrtKpj+vVhlh72lqmmYa/7vd3W+R9E+SFmS7q+idX0maoK45ADslrSizmWya8dckLXL342X20l0PfTXkdSsj7B2Sxna7PUbSwRL66JG7H8wuj0h6Q11vO5rJ4fMz6GaXR0ru5/+5+2F3P+vu5yStUomvXTbN+GuSfuPur2eLS3/teuqrUa9bGWHfKmmimY03s4GSfi5pfQl9fI+ZDc4+OJGZDZb0UzXfVNTrJT2cXX9Y0u9K7OXvNMs03nnTjKvk16706c/dveE/ku5R1yfyn0t6uowecvq6TtKO7Gd32b1Jelldu3Wn1bVH9KikqyRtkvRZdtnSRL39p7qm9t6prmCNKqm3aep6a7hT0vbs556yX7tEXw153fi6LBAE36ADgiDsQBCEHQiCsANBEHYgCMIOBEHYgSD+D1knwahHPak4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# test single image input and check output\n",
    "from matplotlib import pyplot as plt\n",
    "plt.gray()\n",
    "\n",
    "# image number from dataset\n",
    "i = 100\n",
    "\n",
    "inputs = x_test\n",
    "encoded = encoded_imgs\n",
    "decoded = decoded_imgs\n",
    "\n",
    "# inputs\n",
    "print(inputs[i].shape)\n",
    "plt.figure(figsize=(20, 4))\n",
    "plt.imshow(inputs[i].reshape(28, 28), interpolation='nearest')\n",
    "plt.show()\n",
    "\n",
    "# test transpose (for RGB inputs)\n",
    "#plt.figure(figsize=(20, 4))\n",
    "#plt.imshow(inputs[i].T, interpolation='nearest') # TypeError: Invalid shape (3, 128, 128) for image data\n",
    "#plt.show()\n",
    "#plt.imshow(inputs[i].T.reshape(28*1,28).T, interpolation='nearest') # (width * channel, length)\n",
    "#plt.show()\n",
    "\n",
    "# encoded\n",
    "print(encoded[i].shape)\n",
    "plt.figure(figsize=(20, 4))\n",
    "#plt.imshow(encoded[i].T.reshape(8*1, 8).T, interpolation='nearest') # for not flatten\n",
    "plt.imshow(encoded[i].reshape(8*1, 8), interpolation='nearest') # for flatten\n",
    "plt.show()\n",
    "\n",
    "# decoded\n",
    "print(decoded[i].shape)\n",
    "plt.figure(figsize=(20, 4))\n",
    "plt.imshow(decoded[i].reshape(28, 28), interpolation='nearest')\n",
    "plt.show()\n",
    "\n",
    "del inputs, encoded, decoded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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+WH4yJAAATgklEQVR4nO3b+1OU9d/H8Q+6iHIKOaoFtwrlAWvUtNFOZo4w0njKLEfTLE8dNdMcOllRjWlWmmlFOf5ANTSMlifQxgk1NSvSMtM0T4SBCiSgiwjC3n/Anfv9vL6ze/uxeT5+bJ6+ueS6rt1r324hPp/PAAAAAAAAwC2trvQBAAAAAAAA4P9iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCPEkdFRfni4+Ot+9OnT0sH07FjR6mvqamxbv/++29pdrt27aQ+IiJC6hsaGqS2qakpRPoBlxEdHe1LSEiw7svLy6X5Fy9elPq4uDjrtqqqSpodGhoq9cq1bYwxFy5ckPqampoqn89n/8v3w+Px+JS/X3NzszRfvf6V+V6vV5odHh4u9eq9WFtba91eunTJNDc3B+RejI+P93Xu3Nm6P3XqVCB+7GU1NjZat5WVldJs9ZzExMRI/blz56S+rq4uYPdibGysLzk52bovLS2V5qempkr98ePHrduzZ89KsxMTE6W+devWUt+pUyfr9sSJE6aqqiog92KrVq18Ho/941CbNm2k+T6fL2i9+j4UFhYm9VFRUVKvvJ4aY0xTU1PA7sWwsDBfZGSkda++t6ivMwr196a85hijX4PKs1lZWZmprq4OyL0YHh7uU17/1XOo/h7atm1r3R44cECa3b59e6lX/66tWmn/Ll9WVhawe/Gaa67xJSUlWffqZwf1fmlpabFu1fu8R48eUq+8RxujnXev12saGhoCci9GRkb6YmNjrXv1HNbX10u9cj3X1dVJs5X73Bj9GUB1uWdUaWkTHx9vcnJyrPvFixcr483zzz8v9Rs2bLBuP/30U2l2t27dpP7mm2+W+qNHj1q3JSUl0mx/EhISzKJFi6z7V155RZp/6NAhqR8zZox1+9FHH0mzleWUMcZMmzZN6vft2yf1X375pfZpzY/Q0FCjfOA/f/68NL9Xr15Sr7xA7tq1K6jHot6LhYWF1m1FRYU025/OnTtL9/abb74ZsJ/9T5RlwocffijNVs/hqFGjpH779u1SX1RUFLB7MTk52WzatMm6V19nVq9eLfWTJ0+2bvPz86XZ48aNk3p1+fbqq69at/369ZNm++PxeKSFlPLaa4z2jzTGaEvwn3/+WZqtftgfNGiQ1BcVFUl9eXl5wO7FyMhIk5mZad2r7xXbtm2TemVBoDzPGmPMvHnzpP7SpUtSP378eOs2IyNDmu1PTEyM9BrZu3dvaX5TU5PUK58F1GMZMmSI1KvXq7rkmTVrVsDuxaSkJPP+++9b9+oiY+PGjVKvPAMXFxdLs/Py8qReeY82xpg+ffpYt+rvxZ/Y2Fgzd+5c6/7IkSPS/L1790q98o9/mzdvlmanpaVJ/XXXXSf1ytLQGGO+/vrrf7wX+d+jAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzkUeKamhqzZs0a637fvn3SwVRWVkr9xYsXrduPP/5Ymj1t2jSp7969u9S3b9/euvV4pNPkV2Njozlx4oR1v3//fml+r169pP7kyZPW7eDBg6XZxcXFUh8dHS31p0+flvpAiomJMaNHj7buFyxYIM1Xf9edO3e2bsPDw6XZW7ZskfqamhqpnzhxonW7cuVKabY/paWlZsaMGdZ9bm6uNP+bb76R+tdff926HTVqlDT7q6++kvqEhASpV36PxhhTVFQk9f4cO3bMjBs3zrrfvn27NH/ZsmVSr7y+T5kyRZr93nvvSf2tt94q9Zs3b7Zuf//9d2m2P5GRkdKxFhQUSPMHDRok9X369LFur7nmGmn2tm3bpD45OVnqlfclY4xZvny51PsTFRUlvXdNnz5dml9XVyf1ynksLCyUZmdlZUn97Nmzpb53797WbVVVlTTbH5/PJz3b33vvvdL8e+65R+pHjBhh3a5du1aaPXLkSKnv1q2b1GdmZkp9IJWXl5ucnBzrfufOnUE8GmPS09OtW/X31q9fP6l/9tlnpX7dunXWrdfrlWb7c+bMGbNixQrr/tChQ9L8WbNmSb0y/4MPPpBmP/bYY1I/cOBAqQ/U53i+aQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCP+gdCQkKCcRzGGGO2bt0q9atXr7Zuw8PDxaPRHDlyROpLSkqCdCT+eb1e88MPPwRtflZWltQvWrTIum3Xrp16OJIePXpI/a5du4J0JP9ZVFSUGTRokHW/YMECaX5sbKzUL1261Lr1+XzSbPU1Jzs7W+ofeeQRqQ+UkJAQ06ZNm6DNnz59utR7vV7rtqysTD0cyejRo6V+1KhRQTqSK+/8+fNSv3v3but2zJgx6uFI1OtbfQYIlLZt25r09HTrvqCgQJofFRUl9UuWLLFug/16+swzz0j98OHDpT6QamtrTWFhYdDmFxcXS/3Ro0et2/Xr16uHI2lsbJT6ioqKIB2Jfx6PxyQmJgZt/tChQ6V+xowZ1q3ymeS/UVpaKvV9+/YN0pH8Z/X19UH9rBERESH1v/32m3Ub7M8aaWlpUn/o0KEgHYl/HTt2NC+88IJ1P2nSJGl+XFyc1G/atMm6VT8HqFJSUqT+pZdeCsjP5Zs2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOCvH5fPZxSIh9/F8YPHhw0GYfOXJE6keMGCH1y5cvl/qnn37auv3888/N6dOnQ6QfcBnBPocvvfSS1NfX11u3JSUl0uyRI0dK/TPPPCP177zzjjr/J5/P10/6Q5cR7PMYFxcn9W+99ZZ16/V6pdl33XWX1N94441SX1ZWZt1mZWWZffv2XRX34scffyz1S5cutW5nzZolzc7Ozpb66upqqZ89e7bUv/vuu1fNvfjQQw9JfXh4uHXbt29faXZmZqbUp6SkSP2ECROs26KiIlNdXX1V3Ivdu3eX+pkzZ1q3u3fvlmbffffdUj958mSp//DDD6X+0UcfvWruxYULF0p9p06drFv1+WbJkiVSHxKi3SrDhg2zbnfu3Glqa2uvinvxk08+kfoLFy5Yt1VVVdLsV155RerVc1hYWCj1WVlZV829eNNNN0n9Cy+8YN2qz6i//PKL1CvPWsYYM2PGDOt29erVprKy8qq4F+fNmyf1+fn51q16fSi7EGOM2bhxo9TPnz9f6nNycv7xXuSbNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIE8wh48ZM0bq33jjDal/7rnnrNva2lppdkFBgdQPGzZM6k+ePGndNjU1SbMDqbCwUOp37twp9QcOHLBup06dKs3esWOH1D/22GNSn5ubK/VX0n333Sf11dXVUl9VVWXdnjlzRpr9xBNPSP3cuXOlPjs727r966+/pNmBdPDgQalfvHix1Le0tFi3K1eulGa//fbbUh8aGir1+/fvl/orafDgwVIfExMTtH7hwoXS7F9//VXqBwwYIPWtW7e2bkNCQqTZgaT+vc6fPy/1zc3N1u3Fixel2UuWLJH6nJwcqd+8ebPUX0nq61J6errUx8fHW7fqc94nn3wi9V27dpX6u+++27pVXxcC6emnn5Z65ZwYY0xeXp51O2nSJGn27NmzpV69XqdMmSL1V9Ltt98e1Pl79uyxbg8fPizNHj16tNQfP35c6pV73eMJ6sd6vyZOnCj1rVpp3xtJSUmxbi9cuCDNPnv2rNQvWrRI6rds2SL1l8M3bQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHORR4oSEBHP//fdb97m5udLB1NbWSn1aWpp1m52dLc1OTU2V+vz8fKnv1q2bdbtjxw5ptj8ej8fExcVZ93l5edL8jIwMqZ8+fbp1e+rUKWn2ypUrpT4lJUXqb7nlFqn//fffpd6fdu3aSddQSUmJNH/KlClS/+eff1q3TU1N0uyRI0dK/bFjx6R+y5Yt1m1dXZ00258uXbqYN954w7pftWqVNL+4uFjqlXtRPYfq9efxSG9NZvTo0VK/YMECqfcnKSnJTJo0ybpXj/W1116T+qioKOu2Z8+e0uyEhASpj4iIkPqjR49atw0NDdJsfyIiIkyvXr2se/X6nzBhgtQrr0mhoaHS7KFDh0r93r17pb6iokLqA6lDhw7m4Ycftu6vu+46af4NN9wg9RMnTrRuBwwYIM1OTEyU+uuvv17q77zzTutWfdbyp23bttLzd0FBgTRffW/p0qWLdVtfXy/N7t27t9SPGzdO6lVz5swJ2Kzw8HDpNfXbb7+V5s+cOVPqO3XqZN126NBBmq1eg/Pnz5f6a6+91ro9d+6cNNufG264waxYscK6/+ijj6T50dHRUt+jRw/r9rvvvpNmP/DAA1KvfuZfu3at1CcnJ//jf+ebNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADvIocUpKinn//fet+6amJulgWlpapP7OO++0br/77jtp9t69e6V+9uzZUv/4449btzU1NdJsf8LCwkxqaqp1v2zZMmn+jBkzpP6RRx6xbj/77DNpdnFxsdRnZGRIfWJiotQHUlpamlm7dq11/9dff0nz16xZI/UlJSXW7bXXXivN7t69u9TX1dVJ/a233mrdqteUP16v13z//ffWfVRUlDS/TZs2Up+SkmLd7tmzR5rdrVs3qf/pp5+kXn1tCCSv12t+/PFH6z4iIkKaHxsbK/Xt2rWzbpOSkoI22xhjOnbsKPWRkZHW7eHDh6XZ/iQnJ5t3333Xuv/111+l+du3b5d65fofOnSoNNvjkR77TN++faV+/fr1Uh9IkZGRZuDAgda91+uV5kdHR0v9iy++aN22b99emq0+a+3atUvq1depQAkLCzNdu3a17s+ePSvNb2xslPrS0lLrNjc3V5q9YcMGqX/vvfekfuzYsVI/Z84cqfcnKirK3HHHHdZ9dna2NL9///5S/8UXX1i3w4cPl2a3bt1a6seNGyf1ynXS3NwszQ4k9ZmzvLxc6pVn1MrKSmm28gxnjP6MquxO/OGbNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADvIo8f79+01aWpp1n56eLh1MU1OT1B86dMi6XbZsmTRb+XsaY8y3334r9f3797duKyoqpNn+tGrVykRERFj3kydPluavX79e6p988knrdvz48dLs+fPnS/3ChQulPi4uTuo/++wzqffn8OHDJiMjw7qfM2eONF+9nsvKyqzb3bt3S7NPnDgh9f369ZP6X375xbq9ePGiNNuflpYW09DQYN0fPHhQmv/UU09JfW5urnW7detWafasWbOkfs+ePVJ/1113SX0gRUdHmyFDhlj3LS0t0nzl+jTGmAULFli3w4cPl2bn5eVJfc+ePaW+ubnZug0NDZVm+3Ps2DHz4IMPWvdTp06V5qvn8NKlS9btqlWrpNmZmZlSr77mKe8FgdbQ0GD++OMP675Lly7S/M2bN0t9fHy8dTtgwABp9rp166R+x44dUn/hwgXr9rbbbpNm+5OammrWrFlj3S9dulSa37t3b6kPCwuzbteuXSvNVp9VcnJypH7s2LFSH0gej0d6Rq6vr5fmv/POO1KvPIM8/PDD0uyXX35Z6s+cOSP1gwcPDtpsf/7++2+Tn59v3aufG1JSUqReeV9UP2e89tprUq9ef3369JH6y+GbNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADgrx+Xz2cUhIpTGmNHiHg8v4H5/PlxCIQZzDK4rzePXjHP47cB6vfpzDfwfO49WPc/jvwHm8+nEO/x3+8TxKSxsAAAAAAAD8/+B/jwIAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABw0P8CK18jo42jvxgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZQU1RXH8TeJihJQNoEoAooiKKAYUZbIIoKyKCIBNS6oGNGcuMWwKdGoKMclkSMq4oZKjhLQsMRdFmVzBUTRgKAwirLLHkFNJn94cvN7l+lhBrpnqnu+n79u+R41ZVe/quo6776bV1BQEAAAAAAAAJAsPynrAwAAAAAAAMCueGkDAAAAAACQQLy0AQAAAAAASCBe2gAAAAAAACQQL20AAAAAAAASiJc2AAAAAAAACbRPSTrn5eVRH7yMFBQU5KVjP5zDMrW+oKDg4HTsiPNYdhiLOYGxmAMYizmBsZgDGIs5gbGYAxiLOaHQschMG6D05Jf1AQAIITAWgaRgLALJwFgEkqHQschLGwAAAAAAgATipQ0AAAAAAEAC8dIGAAAAAAAggXhpAwAAAAAAkEC8tAEAAAAAAEggXtoAAAAAAAAkEC9tAAAAAAAAEmifsj4AlE+VKlWyuEePHlHbwIEDLa5Zs6bF++23X9Rvzpw5Fp9//vlR2/bt29NynAAAAAAAlBVm2gAAAAAAACQQL20AAAAAAAASiJc2AAAAAAAACcSaNsiYn/70p9H2sccea/FTTz1lccOGDaN+++yzT6FxXl5e1K9Lly4Wv/nmm1FbmzZtLN65c2dJDhtl6KCDDoq2//rXv1r86quvWjxq1Kio37///e/MHlg5petIDRkyxOILL7ww6jd16lSLr7nmmqjt+++/z9DRAQAAlA39XaLx/vvvH/WrWLGixT179oza6tata/Gnn35qsT5XhRDCqlWr9u5gkfWYaQMAAAAAAJBAvLQBAAAAAABIINKjkFZVq1a1eOLEiVFb8+bNLT7ggANS7kPTmXQ6oE+dOfDAAy1u1qxZ1NauXTuLX3vttd0dNsqQTint2rVr1KbndfTo0RaTDpUZPgWxc+fOFg8ePNhiP/W3Tp06Fo8bNy5q86mLyIz69etb3LZt26hNz0F+fn5a/66/LuvY3L59e9RWUFCQ1r+NXekYvuCCC6K2e++91+IxY8ZYPHTo0Kgf19fd89fKWrVqWdyyZUuL/fjQNN/Vq1dn6Oh25Y/jhhtusPiuu+6y2I9ZIFfpGNalGH72s59F/Y4++miLq1WrFrXpMg2nnXaaxXo9CCFONa9QoULK41i4cKHFH3/8cdSP9Cgw0wYAAAAAACCBeGkDAAAAAACQQGWaHuWnlyqmUWenU045xeLjjz8+avvJT/7/jnD27NkWDxgwIOq3ePFii7XyjK6wHkIIM2bMsLhmzZpRm1akIj0q2XSqaL9+/aK2bdu2WUyaTebVqFEj2n7iiScs9ilRat9997V44MCBUds777xj8Y4dO/b2ECG0Qt+1115rcevWraN+7777rsV6393T+6yeb18tbOXKlRb7VLlvv/12j/4eiq9KlSoW33///VGbpi/36dPH4ptvvjnqR3rU7h188MHR9siRIy3Wa+Vnn30W9dP7mD4T/ec//0n3IYZ69epZPH78+Kjtm2++sfi7775L+9/OdnpuDjvssKjtpptuKrTfH//4x6gf6SzJ4n9zVq9e3WIdv926dYv6FbWcQ6rqUV5R91r9d1pld9KkSVE/TcXyqVMoWlHnKRPX3kxhpg0AAAAAAEAC8dIGAAAAAAAggXhpAwAAAAAAkEClvqaN5ltrmcEQQmjQoIHFfm2FWbNmWfzAAw9YrHm5IcRrJhSVQ6hrpfzwww9R257kt5XXNXh8bqCW654zZ07UtnTpUos193fr1q3F+lu+n66r4HOy169fX6x9ouw1adLEYr0GhBDC9OnTLdb1bZA+us7FM888k7JN+eudXgdOPfXUqO2hhx6yWNfO0LVPsGd0XY0TTzzRYl+2V9dWSMe9StehOuGEE6K2vn37WvzKK69Ebaxpk37+Hqzrgh144IEp/92GDRssZg2b4tGywPocGkIILVq0sPjBBx+0eOzYsVG/devWWZzutRR8WW/920ceeWTUptdifR4uz/T89ujRw+K//OUvUb+f//znFutnp2MqhBCGDBlicTatm5Ft9LyFEP/O3Lx5c8p+p59+usU9e/a02JfkVn6s6G8PvY7659X8/HyLp0yZErXpvbtp06YW+9LjX3/9dcrjKq/0/qfXvyuuuCLqp88l/je/Ph+99dZbFvs1+VasWGGx/81ZWu8AmGkDAAAAAACQQLy0AQAAAAAASKBSSY/S6UuVK1e22KdC6BSxihUrRm3a9+yzz7ZY03FCiKfzazlUPy1bp675qXC6T51aV1S/K6+8MmqbP3++xbk8LdJPCXv//fct9tPT9Bz46ftKvy9akk9LB4cQpwb478Enn3xS1GGjDPnp/Dot1ZdW1JK15TUFMR38Z37EEUdYPGrUKIvbtGkT9UtVkrao0pb77bdftH3RRRdZfNppp1ms5z2E+JrJuS6cno8Q4tKkdevWtfiyyy6L+m3ZsiWtx6HXW3++a9asabHe70OgBG4m+OeS/v37W+y/Lzq1/5ZbbrE4l59R0ql3794Wt2/fPmqbPXu2xaNHj7Y43WMvhPj6q+e4bdu2Ub9mzZpZvHz58qjtueeeS/txZRs/PjTVU0tA+/Lu+vnr74zWrVtH/SpVqmSxT+/nHrd3dHkEfa4IIYT69etbPGPGDIt9etGSJUss/uyzzyw+6qijon5r1qyx2KfKzZs3z2K9F/7rX/+K+i1evNhif03Q79P+++8fUinuUhK5RsfY4YcfHrX97ne/s/j888+3WH87ev5+17hxY4v1uzR06NCo38aNGy326bF6zdfvS7ox0wYAAAAAACCBeGkDAAAAAACQQLy0AQAAAAAASKC8kuRV5uXl7XUSpubuadm8EEK49NJLLa5Xr17UpqXUDjnkEIvr1KkT9dPcU80n9TlsmgeuuZEhxKW8NHdfy8iFEOek+ty3e+65p9B+e6qgoCD1IhIlkI5zuJv9FxqHEH8ORX0mmtM5YcIEi7t27Zryby1atChqO+WUUyzWdYnK2LyCgoITd99t9zJ9HjPJ55HrOkiaixxCfB1ISj5vtoxF5dcWmTZtmsXNmze3WHOHPb2G+vGr10xfElOvwzpmfd6vln7/5ptvUh5HmmTlWPQlfSdPnmyxfs6//OUvo347duxI63Fo6dRly5ZFbVqmVNdOCiH9Yzgbx2K6+Weljz76yGL9ToQQwurVqy1u2LChxb48bSlL7FjUNRJDiNet8eXUtUS0rs+VaXocM2fOjNr0mbVXr15Rm67FkQ7ZOBb92pm6/omusemfZbVksP428dauXWvxsGHDorYnnniiWPsoZYkdi/4cdOzY0WJd+zCEuPS6rnPi17TR5xgd6/obM4QQvvjiC4v9tbK4v+9Kcw2jbByLvhx7o0aNLNZ1hFq1ahX109/y+tvCl+TWNYvWrVsXtR133HEW6/OLPyZ9Bvbrsr7wwgsWX3fddRbrd7GECh2LzLQBAAAAAABIIF7aAAAAAAAAJFCplPxWOkXMT1W78847i7UPnSbnpy/p9PFatWpZ7FMydEqxnyr15ZdfWtyyZUuLH3300aiflieeM2dO1FZey/kVNwWqKLVr17ZYSyj6KaQ61V7Ll4awa4oGksOXz9S0GD/lkFK06XHxxRdH2zod1F9DlX7+mirqp9ZrOqifPnzXXXdZrKlY/nug/a644oqorbxeT0OI73dXXnll1KZjR1N0050O5WlZ70MPPTRq0+uyH89ID/1O9O3bN2rT6d3eO++8Y7EvSYtdXXvttdH2YYcdZrGmmoUQwooVK0rjkEIIcWnhO+64w+IGDRpE/V577TWLFy5cmPkDyzJ+iYamTZtarPccLfUbQpzyoKWFfTpd3bp1LfYlgpcuXWrxG2+8UYKjLp98avCtt95qsY7LEEJ49tlnLdZxWtTz5KZNmwqNd/fvsOd0KQx9hgwhhIsuushivadpamII8XuEsWPHWuzHm6bc62/3EELo0qWLxe3atbO4TZs2UT991vH7aNGihcU1atQo9O+GsPfPssy0AQAAAAAASCBe2gAAAAAAACRQqadHFWVPVuH2K0RrqpPGfuVxrbDg/672/eqrryz2U851ZfjSrBaQa3zFmt69e1uslRH8FEVd6X/q1KlR27fffpvOQ0Qa6fTDEOL0HF3hPQTO497QVfWvueaaqE2n1ys/9fTzzz+3+E9/+pPFU6ZMifppqoW/1mpq1j/+8Q+LfdWbPn36WKxV40KIp/mXN40bN7b4xhtvjNr0nvTMM8+U2jFp5ROfepyfn28x08ozQ6dma9XNEOLz4dOENZWGc1M4vX7pNSmEuNqQT63RyofpqH6nz0X+ej1y5EiL9frqq7MNGTLEYn9tR5x+H0L8ma9atcriM888M+qnzylaqW/ixIlRv6Kq1Pbr18/iWbNmWZygSlKJcvXVV0fbmmrtUz1HjRplMd/75NBUwhBCePrppy3u3Llz1KbXYU1P9Cni+ixa3GUx/L3v7bfftljfKfhnm3PPPddiv6yA/lbV1P8lS5YU65iKi5k2AAAAAAAACcRLGwAAAAAAgATipQ0AAAAAAEACJWpNm0wqSZktzRcfM2aMxVWqVIn6jRs3zuJMl1jNNZqvqKXdQohL+WlOoS8l/Mgjj1jsS8uW5xLBSaT5n9ddd13UpjnHw4cPj9pYd2HPaV5tnTp1UvbTz/iDDz6I2rp27Wrx+vXrLS5qfPk2XRfnqaeesnjQoEFRv0qVKlmsJTtDCOGoo46yOB3rRSSZz6PWdRJ8OWctk+nXs0i3WrVqWdy2bVuL/Rj1pTaRfrrOkS+5rtasWRNtL168OGPHlIt8We9GjRpZ7Mfi3//+d4uffPJJi/141tLFWq5Wy9+GEK+Zc8YZZ0RtWkpa+bW/9NqLH+mz5w033BC16bnSNcI+/PDDqJ9e87Rct66DE0K8bpu/L/7iF7+wWNcsYh2//6tdu7bF/lzpM6VfY29vnxF47syMjh07ptz2a8ToOkV//vOfLX799dejfqnWLPJrK+rven/9PP300y3u27evxQ0aNIj6Va5c2WL/Hfniiy8sXrRoUaHHlA7MtAEAAAAAAEggXtoAAAAAAAAkULlJjyqKn5alU+10Wr6WHQshhGHDhllMOk7J6LSz+++/P2rTMolafm3gwIFRv+nTp1vMdMZk0/ScJk2aRG2a9qbnFHvn+OOPt9hf49SKFSssbt++fdTm0w73hJZhfO655yzu379/1E/TDXwqardu3SweO3bsXh9TkmkaUgjxFHFNUQshnjac7nuQT+vo1KmTxVoad9OmTVG/SZMmpfU48COd7q3lUfVchBB/D/xUcl8aF7vSz++yyy6L2ubOnWuxpp+GEE+5v+WWW1LuX0s6ayrMzp07o36aMqOpoyHE3wUdfwMGDIj6Ue54V/rZ+c9Vr3k+NSLVPqpWrZpyf6n+TWHb2NXMmTMt1pLKIcTfbU27DmHP7oX8hsi8mjVrRtt67/LnbN9997W4V69eFh9yyCFRv7Vr11qsY7Zly5ZRvxo1aljsU1H1b2lc1Bj199IRI0ZY7J+J0omZNgAAAAAAAAnESxsAAAAAAIAEIj0qxGkEIYTQunVri3XK1iWXXBL1y/UqJulWsWJFi59//nmLfRUGnaao05G1IkMITP1NOp1aqBWjdAX2EEKYP3++xVu2bMn8gZUTxx57rMU+PUrHjk49TUc6VFE0FU6ntYZQdBUcP602l/lp+VqZcM6cOVGbT6lIJ01TDSFOUdMUj/fffz/qt3nz5owdU3mm6TI6Zn0am6YjPv7441Ebadwlo6mjIcRjU6uMhBDCr3/9a4t1Cr+vBqQVnebNm2fxhg0bon79+vWzuGnTpimP8frrr7d45cqVKfvhR0VVS6xfv77FqSrKhJC6GqZWBvN/y9PqYJpqp1VoyiO9nlWvXt3iolJV/PjQZ0pNYynu9a8kqVJ6XFxfi/a3v/0t2u7evbvFJ510UtSm1Z50SYVmzZpF/fz973982nBx0xH1HOpzTggh5OfnW3zhhRdGbe+++26x9r+3mGkDAAAAAACQQLy0AQAAAAAASCBe2gAAAAAAACRQuV3TRvPgzjvvvJRtM2bMsPjll1/O/IHlEJ9r2Lt3b4u1lLrPG9Sytloi2PdDsmlO6TnnnGOxz/u99957LeYcp88xxxxjsc/n1TzvZcuWZfQ49G8fffTRFmu+uue/B1OnTk3/gSXU8uXLo21d68Kv7aPrgW3durXEf8uvdaTrjrVq1Spq69ixo8VaFnPcuHFRP8ZwZmjp96LWONH7p66ZEgJrLuwtXZ/m4Ycfjtr8dkn5NRjq1KljcePGjaO2jz76yOLx48fv1d8tzy6//PJo++STT7ZYr7X+3Oq50vVP/Hj78ssvLe7QoUPUpiWs77jjDov92pnl7Xqq9yR9TtHS6iHE58CXutfz+PHHH1u8cePGqJ+uoahloPV3Xwjx2nH16tVLuY/PPvus0H+DH61bty7aPvPMMy3297SzzjrLYl13tlq1alE/vd/p2lC6PlUI8f3TP/foGNb96XpVIYQwefJki33J79LCTBsAAAAAAIAE4qUNAAAAAABAApXb9CidVq5pOyHEJTN1elRJysAhhFq1akXb9913n8VavvTrr7+O+ulUR02P4vPPLlr+UlNh/BTJ6dOnl9oxlScnnHCCxT49Sqf5p3tc+bRIPY4HHnjAYl/iW49jyZIlUdunn36azkNMtG+++Sba1lKSWnY7hLhkrZYDf+mll6J+mjqlJTP9uWrRooXFPiWjRo0aFuuUfVJwMsOP2VNPPdViTU/zn/e0adMsZop+9vBpMOeee67F/ho9aNAgi8tqmn4u2LRpU7StqUmadla5cuWU+5g7d67FXbp0idr0vA0ePDhqu+mmmywuKk3ElyXPdfr7S9ONfDp1hQoVLD7kkEOiNk3H79Gjh8Xfffddyn0on0alqXI+tUbH7fDhwy2+7bbbon78ftmVnmst0x5CCAsWLLC4qLLq2qa/OR588MGoX58+fSz250LPt6ZMvvjii1G/JJxDZtoAAAAAAAAkEC9tAAAAAAAAEoiXNgAAAAAAAAlUbta08fnhWmbYr63w/vvvW7x48eLMHliO0c/5hhtuiNo03/CHH36w+Pbbb4/6aS4x6yNkr7Zt2xb63/Pz86NtLbGHPeevcT7PW+maGP7f7Qktlzls2LCorV+/fhbr2gD+727bts3i3//+91Hbjh079voYs4WuNxRCnCfv1zs44ogjLO7Zs6fFvXr1ivppLr9+7v5z1fKlVapUSXmMev324xnp4cfHxRdfXGg/vxbKyJEjLeb+mT20rG0Icclav8bG22+/XRqHVO7MnDnTYl3P8le/+lXUb9y4cRZPmDDBYr0ueo899li0rWtn6Bqb/tq9aNGiYu0/V+g1q3v37oXGIYTQt29fi/3aNPrso7/vKlWqFPXTa6yu7+Z/E/q131K16XOLrnUUQgivv/66xVyXd08/o+J+XnoN1eehEOIS8f6eOWbMGIt1PcAkrGHjMdMGAAAAAAAggXhpAwAAAAAAkEB5JZmmlZeXl7Vzutq1axdta5lhPwWqefPmFuvUxLJUUFCw9zkMIfPnUKcbvvfee1Gblv5dtWqVxX7Kvy95W5y/lenphlqi3P8tnWq3m+l08woKCk5Mx/EkcSz66fzPPvusxTrl9+GHH476XX311Zk9sDRL6lj0n7+mnVWrVi1q03QZLQFdVGltnV7apEmTqO3xxx8vdH/+3+kx+imq+n3p379/1JaBsrZZORb1OhRCCHXr1rVYp8771CZ/Tv5n3bp10fby5cst9mluZ599tsWawlW7du2o3/bt2wv9W5mQ1LGYDpoyEUIIK1assFjL365cuTLqd+SRR1rsS9wmVFaOxXTQ6+HSpUujtgYNGlj87rvvRm1t2rSxOCkpM7k2FotKidmTtAm9D4YQwpAhQyzW+92GDRuifp06dbLYX68zIGvGoo4dTfcOIU4tHDhwoMUdOnSI+h144IEWF5W6rc8q/rlFU8OV//3Tvn17izN9Xc61saj8uWnVqpXFmoJWsWLFlPvQUvIhxL9BfXp6GSp0LDLTBgAAAAAAIIF4aQMAAAAAAJBAOV09SlcQf/XVV6M2nfq4cOHCqO2TTz7J7IGVE36avE7j1emlfur+ggULLNZUJD9dXKcp+ioomk6h+/DTKDXd4JhjjonaLrjgAot79Ohh8QEHHBD1+8Mf/mDxnDlzorZly5aF8sJP/9XqUTqlUVdnR+asWbPGYp8epde/Rx991OJnnnkm6qdVFLRaQ61ataJ+++zz/1uJn1au40/HrE/ruOWWWyzOQDpUTvDTqot7ffnggw9K/LdmzJgRbZ9zzjkW6zn14x7pcdxxx0XbWn1R75+33XZb1C9LUqIQ4ucRf03V8/jII49EbT5FA+mX7soxfn96r9UUKP89uOqqqyy+++67o7byVFXR03uQv+Zpmvf1119v8YgRI6J+Z5xxRqH73rlzZ7Stvwm///77qE1/N/jfKIr7ZHr4dOwnn3zS4qJSorQ6qaZ6h5ColKjdYqYNAAAAAABAAvHSBgAAAAAAIIF4aQMAAAAAAJBAObemjeYNTpkyxWJfKlXLSvu8xnTnspYnmmf68ccfR21amq1GjRoWP//881E/PVcVKlSw2Jd60/Pkc1q1VLuuQVO1atWon7b5tWq0lF+qssUhhDBy5EiLtfRxCHE+ba7zawLVrFnTYv1e/POf/yy1YypPfCn6Sy+91OI333wzatMxpuOydevWUT9dn8Z/71Px10/ND9drwrXXXhv1y8/PL9b+UTq0HGoIqcu1s75GZlx33XXRtq4bpWPqlVdeKbVjQnqddNJJFvu1wLQs7bRp04q1v6Ku0f7+gKLpZ5mOz87vY/ny5RbrujUPPfRQ1E/XVpw5c2bU9sYbb+z1ceW6rVu3WjxmzJio7ayzzrJY15fyY7FRo0YW+/VP9Lqs69289dZbUT+/Tg6KT59X/bpOhx9+eKH/xv8m1Pup/j7MNsy0AQAAAAAASCBe2gAAAAAAACRQ1qdH+emgmuqk6RpbtmyJ+nXo0MHitWvXZujoyrcNGzZE2zo9VNOeNA6h+GkYypfybtGiRaF/10/l1zadRhlCCF9//bXFlStXttiX9XvwwQctXrVqVUkOO6fcc8890bamlGl5Z1/qGZmhZZ5ffPHFqE2nBev03uLyU71/+OEHizdt2hS1LViwwOKbb77Z4nnz5kX9SLPJTkz7zgxf8lvvi3qvWrduXakdE9Jr+PDhFvuUDH3+8M+v2levxX4fmqpKelTJZPrz0v1rierJkydH/a655hqLe/fuHbXNmjXLYu6fu+efOTQFUdOBfXnuSpUqWazLJoQQjzFN/R81alTKfiiZrl27Wty9e/eoTc+VjgGfRvXEE09k6OhKFzNtAAAAAAAAEoiXNgAAAAAAAAnESxsAAAAAAIAEyvo1bapXrx5t+zzC/5k6dWq0/cknn2TsmPCjRx99NNr+/PPPLdYSh7Vr1476HXTQQYXuz+fsakk3n9evZfl0jY358+dH/V599VWLly5dGrXpv9N1G7SsX3mnpfh8yW/N4R0xYoTF5F6XDv3ODho0KGrT8vYdO3a02K8NlarsqV/HRNdymjRpUtT29NNPW6wlv3UdHCRPp06dUrZpeXbGc/roeNuxY0fKfrouGGslZBe9xh577LEW+7X9mjRpYrFfy+SFF16wWJ+DfJnbbdu27d3BolTovVDvlyGEMGDAAItPP/30qO2oo46yeMmSJRazflHh/DX1yiuvtPixxx6zuGnTplE/XSvKl/zW8tGXXnqpxbpeDkpOr4e33367xbr2kPfhhx9afOedd0ZtuTImmGkDAAAAAACQQLy0AQAAAAAASKCsTI/S6aU+HerQQw+1WKcc3nfffVE/pnRn3hdffBFtjxkzptAY2UvH4ldffRW1aclSnXqaK9MUs8myZcui7csvv9xiTY/S6b0hhHD00UdbrGWGx48fH/V76aWXLNYUqBBC2L59u8WkciSbls88+OCDozZNC3355ZdL7ZjKE02PWr16ddTWuHFji3VM+fK0pO8mW8WKFS3eZ5/Uj+BVqlSx+IILLoja3nvvPYs3bNhgsU+x0u8JssPGjRuj7blz51p80kknRW2a9vzb3/7WYp/Cgx/5Z09Np7nkkkssvvHGG6N+Oo6efPLJqG3BggUW6+fOc+7eGTp0qMWNGjWyWFPVQog/cz1vuToGmGkDAAAAAACQQLy0AQAAAAAASKC8kkzhysvLS8R8r86dO1vsp2nr1CmtKNS2bduo3+LFizN0dJlRUFCQt/teu5eUc1hOzSsoKDgxHTviPJYdxmJOYCwWQtMdx40bF7Vplb/zzjvP4i+//DLzB5ZCLo/Fnj17Rtt33323xYMHD7Z44sSJUb8sTEEsV2NRU+BatWplsU/hX79+vcWjR4+O2l577TWLi0r11yUCMp2ukctjsSxVq1bNYl/1Vqs2aspc//79o34lOPflaizmqmwciz5VVL/rWiXNf5e1MnT37t0t9pX0slChY5GZNgAAAAAAAAnESxsAAAAAAIAE4qUNAAAAAABAAmVlye8ePXpYrPnB3ooVKyxeu3ZtJg8JAICspuWie/XqVYZHAr9Wjd9GdtI1GbSc88knn1wWh4OE27x5s8X62yeEEBo2bGjxhAkTLKbcNLJNlSpVUrbpul26Vm0IIVx//fUW58A6NrvFTBsAAAAAAIAE4qUNAAAAAABAAmVletTs2bMt/s1vfhO1bdu2zeJu3bpZvGnTpswfGAAAAADsJU0Neeedd6K2rVu3Wrxjx45SOyYg3TQNMIS4bP1VV11l8a233hr109Lg5QEzbQAAAAAAABKIlzYAAAAAAAAJxEsbAAAAAACABMorSWm4vLy8xNWRK6rkdy6VvSsoKEj9P1oCSTyH5ci8goKCE9OxI85j2WEs5gTGYg5gLOYExjH8B5YAAADGSURBVGIOYCzmBMZiDmAs5oRCxyIzbQAAAAAAABKIlzYAAAAAAAAJVNKS3+tDCPmZOJA9lUspUEWol8Z9Je4cliOcx+zHOcwNnMfsxznMDZzH7Mc5zA2cx+zHOcwNhZ7HEq1pAwAAAAAAgNJBehQAAAAAAEAC8dIGAAAAAAAggXhpAwAAAAAAkEC8tAEAAAAAAEggXtoAAAAAAAAkEC9tAAAAAAAAEoiXNgAAAAAAAAnESxsAAAAAAIAE4qUNAAAAAABAAv0X9Dff6/bupMMAAAAASUVORK5CYII=\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.99 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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0014, PSNR (dB) = 28.68\n",
      "MSE = 0.0049, PSNR (dB) = 23.07\n",
      "MSE = 0.0007, PSNR (dB) = 31.69\n",
      "MSE = 0.0030, PSNR (dB) = 25.25\n",
      "MSE = 0.0026, PSNR (dB) = 25.92\n",
      "MSE = 0.0006, PSNR (dB) = 31.88\n",
      "MSE = 0.0048, PSNR (dB) = 23.17\n",
      "MSE = 0.0036, PSNR (dB) = 24.49\n",
      "MSE = 0.0048, PSNR (dB) = 23.21\n",
      "MSE = 0.0035, PSNR (dB) = 24.56\n",
      "Wall time: 0 ns\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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9834\n",
      "0.9632\n",
      "0.9903\n",
      "0.9734\n",
      "0.9687\n",
      "0.9924\n",
      "0.9357\n",
      "0.9684\n",
      "0.9584\n",
      "0.9698\n",
      "Wall time: 200 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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9983\n",
      "0.9889\n",
      "0.9967\n",
      "0.9963\n",
      "0.9944\n",
      "0.9988\n",
      "0.9845\n",
      "0.9922\n",
      "0.9918\n",
      "0.9916\n",
      "Wall time: 744 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",
    "    d=MultiScaleSSIM(original, contrast, max_val=1.0)\n",
    "    print(f'{d:.4f}')"
   ]
  },
  {
   "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": 4
}