Simple Autoencoders in Keras - 784-64-784-tf1-rpi.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip3 install -U matplotlib opencv-python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoding_dim = 64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_img = Input(shape=(784,))\n",
    "encoded = Dense(encoding_dim, activation='relu')(input_img) # apply to previous layer = input_img\n",
    "decoded = Dense(784, activation='sigmoid')(encoded)\n",
    "\n",
    "autoencoder = Model(input_img, decoded)\n",
    "autoencoder.run_eagerly = True # plaidml / keras not running eager mode by default\n",
    "\n",
    "encoder = Model(input_img, encoded)\n",
    "\n",
    "encoded_input = Input(shape=(encoding_dim,))\n",
    "\n",
    "decoder_layer = autoencoder.layers[-1] # extract last layer\n",
    "\n",
    "decoder = Model(encoded_input, decoder_layer(encoded_input))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.compile(optimizer='adam', loss='binary_crossentropy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1595257695\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 43s - loss: 0.6924WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.160000). Check your callbacks.\n",
      "235/235 [==============================] - 18s 79ms/step - loss: 0.2440 - val_loss: 0.1595\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 19s 81ms/step - loss: 0.1415 - val_loss: 0.1249\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 16s 70ms/step - loss: 0.1162 - val_loss: 0.1066\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 18s 77ms/step - loss: 0.1020 - val_loss: 0.0959\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 15s 63ms/step - loss: 0.0933 - val_loss: 0.0890\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 17s 71ms/step - loss: 0.0875 - val_loss: 0.0843\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 21s 88ms/step - loss: 0.0835 - val_loss: 0.0811\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 17s 72ms/step - loss: 0.0808 - val_loss: 0.0789\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 17s 71ms/step - loss: 0.0789 - val_loss: 0.0774\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 17s 71ms/step - loss: 0.0776 - val_loss: 0.0762\n",
      "Epoch 11/50\n",
      "235/235 [==============================] - 15s 64ms/step - loss: 0.0766 - val_loss: 0.0754\n",
      "Epoch 12/50\n",
      "235/235 [==============================] - 30s 129ms/step - loss: 0.0759 - val_loss: 0.0749\n",
      "Epoch 13/50\n",
      "235/235 [==============================] - 43s 181ms/step - loss: 0.0753 - val_loss: 0.0744\n",
      "Epoch 14/50\n",
      "235/235 [==============================] - 27s 114ms/step - loss: 0.0749 - val_loss: 0.0741\n",
      "Epoch 15/50\n",
      "235/235 [==============================] - 27s 116ms/step - loss: 0.0746 - val_loss: 0.0739\n",
      "Epoch 16/50\n",
      "235/235 [==============================] - 28s 121ms/step - loss: 0.0744 - val_loss: 0.0736\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 27s 115ms/step - loss: 0.0742 - val_loss: 0.0735\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 28s 120ms/step - loss: 0.0740 - val_loss: 0.0734\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 30s 126ms/step - loss: 0.0739 - val_loss: 0.0732\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 26s 112ms/step - loss: 0.0738 - val_loss: 0.0732\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 30s 128ms/step - loss: 0.0737 - val_loss: 0.0731\n",
      "Epoch 22/50\n",
      "235/235 [==============================] - 30s 127ms/step - loss: 0.0736 - val_loss: 0.0730\n",
      "Epoch 23/50\n",
      "235/235 [==============================] - 31s 133ms/step - loss: 0.0735 - val_loss: 0.0729\n",
      "Epoch 24/50\n",
      "235/235 [==============================] - 32s 137ms/step - loss: 0.0734 - val_loss: 0.0728\n",
      "Epoch 25/50\n",
      "235/235 [==============================] - 31s 134ms/step - loss: 0.0734 - val_loss: 0.0728\n",
      "Epoch 26/50\n",
      "235/235 [==============================] - 30s 129ms/step - loss: 0.0733 - val_loss: 0.0727\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 32s 137ms/step - loss: 0.0733 - val_loss: 0.0727\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 33s 141ms/step - loss: 0.0732 - val_loss: 0.0727\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 36s 154ms/step - loss: 0.0732 - val_loss: 0.0727\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 35s 150ms/step - loss: 0.0732 - val_loss: 0.0726\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 36s 153ms/step - loss: 0.0731 - val_loss: 0.0726\n",
      "Epoch 32/50\n",
      "235/235 [==============================] - 35s 148ms/step - loss: 0.0731 - val_loss: 0.0726\n",
      "Epoch 33/50\n",
      "235/235 [==============================] - 38s 163ms/step - loss: 0.0731 - val_loss: 0.0725\n",
      "Epoch 34/50\n",
      "235/235 [==============================] - 39s 164ms/step - loss: 0.0730 - val_loss: 0.0726\n",
      "Epoch 35/50\n",
      "235/235 [==============================] - 39s 166ms/step - loss: 0.0730 - val_loss: 0.0726\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 38s 162ms/step - loss: 0.0730 - val_loss: 0.0725\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 38s 160ms/step - loss: 0.0730 - val_loss: 0.0725\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 36s 153ms/step - loss: 0.0729 - val_loss: 0.0725\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 36s 151ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 37s 157ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 32s 138ms/step - loss: 0.0729 - val_loss: 0.0725\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 36s 154ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 38s 162ms/step - loss: 0.0729 - val_loss: 0.0725\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 35s 148ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 39s 165ms/step - loss: 0.0728 - val_loss: 0.0724\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 37s 156ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 37s 158ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 34s 146ms/step - loss: 0.0728 - val_loss: 0.0724\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 32s 137ms/step - loss: 0.0728 - val_loss: 0.0724\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 32s 138ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "1595257695\n",
      "Wall time: 25min 6s\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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoded_imgs = encoder.predict(x_test)\n",
    "decoded_imgs = decoder.predict(encoded_imgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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+WH4yJAAAUJElEQVR4nO3d+W+UZbiH8Wfa6Uad2kILBVqtgFQqi2GJohKiBQNoRRHEDdmqIFsUVFRiFBcwJYGogRDEyCKoCIqAklBAKCWAVBAUkaW1lQq1QOkytHSd8wccGZ/vyXh8NNfnR3NxO8w77zJ3a/QEAgEDAAAAAAAAt4T90y8AAAAAAAAA/xtLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB3ml2OsNREVFWffq/07c65VejvH7/dZtfHy8NLuurk7qIyMjpb6hocG6bWxsNE1NTR7pX3AVMTExgbi4OOu+urpamq8ew3bt2lm3paWl0uympiapDw8Pl/oePXpI/ffff38hEAgkSX/oKrxeb0D5zLVt21aaX19fL/VJSfZ/raKiImm2eh1pbGyU+uuuu866LS8vN9XV1SE7F30+n3VfVVUlzVePuXLuVlRUSLM9Hu0tS05Olnr181pcXByyczEuLi6gfP7V90K5zxmjfU4SEhKk2bGxsVKv3heVa3ZZWZmpqqoKybkYHh4eUD7/f/f1VLkm/fLLL9Lsv/NZxRj9M1JeXh6yczE2NjagPOupzzfqudilSxfrtqamRprd3Nws9ep9UXlOrKioMH6/PyTnYlhYWCAszP7nyampqdJ89Zir77NCPVfU15KSkiL1oXxGjY6Olp5vWlpapPnqdUl5pjh37pw0+9prr5V69buJ8gxcU1Nj6urqQvaMqvzd1GuMSjnmLl0f/y/za2pq/vRclL5hR0VFmYyMDOtePamUB19jjMnLy7Nu77nnHmn2Tz/9JPXKQ5Yxxvz222/WbWFhoTQ7mLi4OPP4449b99u2bZPmJyYmSv2MGTOs2xdffFGarX6xVG4wxhhTUFAg9R6Pp0T6A0FERkaa9PR0637KlCnS/JIS7aVmZ2dbt0888YQ0W73Yqcu9hQsXWrczZ86UZgfj8/nMyJEjrfutW7dK86dOnSr1yrm7bt06aba6zH3ppZek/uTJk1I/fvz4kJ2LSUlJZt68edZ9dHS0ND8/P1/qt2zZYt2OHj1amt23b1+pV79QXbx40bp95plnpNnBeL1e06FDB+t+2rRp0nz1Hr5kyRLr9s4775RmK39PY7RnFWOMue2226T+3XffDdm5GB8fbyZPnmzd79y5U5q/a9cuqV+0aJF1qzzPGmNMZWWl1P/xxx9SP3jwYOs2JydHmh1MWFiY9EPW1157TZq/Y8cOqVfeZ/UHf3369JF6dbGnHpdQPqP6fD7z0EMPWffqMu3s2bNSP2vWLOv2nXfekWbfe++9Ul9WVib1ypJn/fr10uxgrr32Wul5vby8XJqvLq9+//1361b9Iad6fRw0aJDUq8d8+/btf3ou8p9HAQAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADjIq8S1tbWmoKDAun/00UelF1NRUSH1fr/fuh0+fLg0W9WqVSup79Spk3VbWlqqvpyrqq2tNQcPHrTuu3XrJs1/+eWXpX7s2LHWbefOnaXZr776qtRv3bpV6pOTk6U+lDIyMqRzccaMGdL8lpYWqf/ss8+s28mTJ0uzN2zYIPXquf7FF19Yt5WVldLsYCoqKszatWut+2HDhknzFy9eLPVr1qyxbqOjo6XZly5dkvo33nhD6uPi4qQ+lIqKiqR7XSAQkOZ/8MEHUq+8luLiYmn29u3bpV65RhljzNy5c6U+VG688Uazbt066z4nJ0ea3759e6lXrkmZmZnSbOX+b4wx/fv3l/pQXiNVFy9elK5jQ4cOleaPHz9e6pXnp82bN0uz1WfOm2++WerPnDlj3TY0NEizg4mIiJDOl2+//Vaav2fPHqnPysqybuvq6qTZ8+bNk/oHH3xQ6tVn2lCqr683J0+etO47dOggzVevM/n5+dbt5cuXpdmnT5+W+lOnTkn9W2+9Zd3u2LFDmv1XwsPDrdvExERptvr8sXv3buv2sccek2bffvvtUq+8L8bo14ar4TdtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAc5FXiqKgok5aWZt3HxcVJL2bHjh1SP3DgQOu2X79+0uyjR49KfUFBgdQnJydbt01NTdLsYMLCwkyrVq2se7/fL81fvny51F+5csW6nThxojT7+uuvl/olS5ZI/YwZM6R+zpw5Uh9MSUmJmTRpknW/fv16af4TTzwh9QcOHLBuvV7psiOfi+3atZP6YcOGWbd79uyRZgcTGxtrbr/9duv+lVdekea//vrrUj9lyhTr9pZbbvlbX8vOnTulPhAISP3nn38u9cFERESYpKQk6149t1q3bi31ymc0PT1dmj1hwgSpr6qqkvq8vDzrtqamRpodzB9//GEWLVpk3X/55ZfS/OHDh0v9p59+at327dtXmn348GGpV54XjDEmJSVF6kMpMjLSpKamWvdHjhyR5peVlUl9TEyMdVtYWCjNHjx4sNSXlpZKvc/ns249Ho80Oxj1GJ47d06ar34+V65cad326tVLmj158mSpz8rKkvrMzEypD6X6+npTUlJi3avXyMbGRqn/8ccfrdt77rlHmp2TkyP1s2fPlvrc3Fzrtrq6WpodTFNTkzl//rx1r14fIyIipP7pp5+2btVnyKeeekrqle+uxhhz4403Sv3V8Js2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAO8ipxVFSUSUtLs+4TEhKkFxMWpu2Q2rRpY90ePHhQmv3www9L/ebNm6W+d+/e1m1RUZE0O5S2bdsm9SNGjJD6sWPHWrerV6+WZldUVEj9hAkTpH7//v1SH0p+v9/k5eVZ9xMnTpTmL1iwQOrffvtt6zY3N1eaXVNTI/X5+flS37p1a+u2trZWmh1MU1OTuXjxonU/depUab7y+TDGmDlz5li36jHMycmR+tjYWKnfvXu31IdS27ZtzbRp06z7AwcOSPPV/v7777duw8PDpdnqdeTee++V+pKSEuu2oaFBmh1MfX29KSwstO4feeQRaf6yZcukPjs727otLS2VZkdGRkr9kSNHpD4mJkbqQ8nr9UrPnSNHjpTmL1++XOrfe+8967ZHjx7SbPWZs2fPnlJfWVlp3TY3N0uzg/H5fGbAgAHWvXrt37dvn9T379/fulWfCU+cOCH16rk+Y8YMqVc+r3/F5/OZgQMHWve//vqrNF+5XhtjzKRJk6zb48ePS7OTkpKkfsyYMVJfVVVl3TY1NUmz/2qW8j2qU6dO0vyvv/5a6ouLi63bUaNGSbO3bt0q9RkZGVJ/6dIlqb8aftMGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzkVWK/32/y8/Ot+x07dkgvpk+fPlIfFxdn3Z46dUqaffbsWalvbm6W+suXL1u3LS0t0uxg2rZta6ZPn27dd+3aVZr/22+/Sf2GDRus24aGBml2v379pL6srEzqV6xYIfWh5PF4THR0tHW/detWaX5aWprUR0ZGWrcxMTHSbPUz2L59e6nfv3+/dauct3+lY8eOZu7cudb92rVrpflDhgyR+qVLl1q3SUlJ0uzDhw9L/X333Sf1JSUlUh9KFy5cMB999JF1/+GHH0rzd+3aJfUzZ860bhcsWCDN/uqrr6R+2bJlUq88A6jPF8EEAgETCASs+/Xr10vzk5OTpT4hIcG6LSwslGbX1tZK/dChQ6VevU+HUktLi6mrq7Pu1dc6ZcoUqVeeb9Rr5OzZs6X+oYcekvrs7GzrVjl3/orf75fuyYmJidL84cOHS71y/b3tttuk2REREVKv3GeMMeb555+X+lAKDw+XvqMdOXJEmq/MNsaYgoIC67aoqEiavWjRIqnfvn271Pt8PuvW4/FIs4NpaWkx1dXV1v2+ffuk+W3btpX69PR061b9nqF+h1ff544dO0r91fCbNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADvIqcWRkpLnhhhus+4yMDOnFnD17Vur9fr91e8cdd0iz27VrJ/X19fVSn5eXZ93W1NRIs4MpLS01L7zwgnU/YMAAaf6VK1ekPiEhwbodMmSINPvChQtSrx7D9PR0qT9x4oTUB5OQkGBGjRpl3YeFafvZxYsXS/2ePXus24iICGn2yJEjpX7fvn1S/+uvv1q36mckmIqKCrNu3Trr/ty5c9L8pKQkqe/Tp491q17by8rKpH7ZsmVS39zcLPWhdM0115j+/ftb96tXr5bmnz9/Xur79etn3S5dulSanZ+fL/Wff/651Hfp0sW6ra6ulmYHEx8fb7Kysqz7wYMHS/NXrVol9SdPnrRu09LSpNnqa1evOxs3bpT6UIqMjDQpKSnW/d69e6X5qampUq/ci06fPi3NVu+LTz/9tNT/U9fUxsZG6buA+hzWsWNHqa+srLRun332WWn2uHHjpD43N1fq27dvL/Wh1NTUZMrLy637uro6aX58fLzUHzt2zLp97rnnpNlVVVVSf/jwYan3eDzWrfK9+K9ERERI54v6bK/cc40xZsWKFdbtypUrpdnZ2dlSr34vV77zB8Nv2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOMirxB6Px4SHh1v3zz//vPRiVq5cKfV79uyxbuvr66XZW7ZskfpevXpJ/ezZs63b+fPnS7OD8fl8JjMz07o/deqUND8QCEh9bW2tdXvw4EFpdlpamtS3adNG6i9fviz1oXTlyhXzyy+/WPedO3eW5s+aNUvqv/vuO+v2ypUr0myfzyf1ZWVlUn/XXXdZt5s2bZJmBxMdHW1uuukm615939TPc6dOnazbrl27SrMPHTok9dOnT5f6/fv3S31ubq7UB+P1ek1iYqJ139jYKM2PiIiQ+rlz51q3I0aMkGar9wO/3y/1ynVEfV+CqampMTt37rTue/fuLc1/9tlnpT4/P9+6raurk2anp6dL/ZIlS6R+4MCBUr9u3TqpD6aurs789NNP1v24ceOk+dnZ2VK/cOFC6/b111+XZvft21fqjx8/LvVnzpyxbtXn62ASExPNxIkTrftWrVpJ8998802pj4qKsm7XrFkjzR4zZozUHzt2TOo3btwo9aFUUVFhPvnkE+t+6NCh0nz1+aZ79+7W7bZt26TZgwYNkvrTp09Lfc+ePa3bsLDQ/S5GQ0ODdB1Q7y3Kc5MxxnTr1s26Vb/DqNewpqYmqX///felPisr60//Ob9pAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIK8SZ2RkmIKCAus+KytLejHV1dVSP3LkSOt248aN0uyqqiqpb926tdTPnz/fui0rK5NmBxMfH28eeOAB637atGnS/FtvvVXqz58/b90++eST0uzdu3dL/dq1a6W+W7duUl9aWir1wcTExJgePXpY9/n5+dL8QCAg9e3atbNujx8/Ls2+cOGC1JeXl0t9x44drVv1fQkmOTnZvPjii9a9cr0zxphdu3ZJfWNjo3WrXqu9XulWY0aPHi31nTt3lvrc3FypDyYuLs4MGTLEut+wYYM0PyEhQeq3bNli3c6ePVuarZznxhgzb948qa+trbVuW1papNnBtGrVyvTu3du6V6/lmzZtknrlHv3zzz9Ls4uLi6Xe4/FIvXIdCbW0tDSzatUq637JkiXS/KeeekrqBw8ebN0eO3ZMmv3NN99IfXJystQrx/3o0aPS7GD8fr/Zu3evda9eT7/99lupnzp1qnV73333SbOnTJki9aNGjZL6ESNGSP0PP/wg9cHExsaa7t27W/fqM4Iy2xhjioqKrNtDhw5Js9u0aSP1d999t9SnpKRYt4WFhdLsYJKSkkx2drZ1v3r1amn+xx9/LPWpqanWbWJiojT7k08+kfrMzEypX79+vdRfDb9pAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIE8gELCPPZ7zxpiSv+/l4CquDwQCSaEYxDH8R3Ec//04hv8NHMd/P47hfwPH8d+PY/jfwHH89+MY/jf86XGUljYAAAAAAAD4/8F/HgUAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgoP8BaRKinOB2A1IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZQU1dnH8cubSEARQVQ2RcEFEEFFFHFHVKKIG25gjB5RXON6ErcYxCURjIL7EkHwxCUkQU1UJG64QEQEJeCKIG6gKCKiKGoy7x8eH3/3Yaqo7plhqru/n7+ePvdOdU/dvtXVfe5znwZVVVUBAAAAAAAA+fJ/9f0CAAAAAAAAsCp+tAEAAAAAAMghfrQBAAAAAADIIX60AQAAAAAAyCF+tAEAAAAAAMghfrQBAAAAAADIoZ8W0rlBgwbUB68nVVVVDWrjOIxhvfqkqqpqw9o4EONYf5iLZYG5WAaYi2WBuVgGmItlgblYBpiLZaHauchKG2DNeae+XwCAEAJzEcgL5iKQD8xFIB+qnYv8aAMAAAAAAJBD/GgDAAAAAACQQ/xoAwAAAAAAkEP8aAMAAAAAAJBDBVWPAgAAAACgkv3f//249uF///tfPb4SVAJW2gAAAAAAAOQQP9oAAAAAAADkEOlRWGN++tMf327bbrutxeedd17U78ADD7S4UaNGFv/kJz+J+v33v/+1+M4774zaTj31VItZsli6GjRoUG1cVVUV9fOPAdSMzrcQmGMAACi+X2BNYqUNAAAAAABADvGjDQAAAAAAQA7xow0AAAAAAEAOsacNasTve6Dl7zbeeOOobcKECRZ369bNYr9XTdrxk9r23HPPqG2dddaxePny5YnHKIZ/TbpXz7ffflurz1UJ9PwdcMABUZvub3TPPfdYPGXKlKgf571u6J5S11xzjcVHH3101O+jjz6y+IQTTojapk2bZnHWfVHYT6Vm0q6bWc8l57z06edxCPG19rvvvrOYfRkAYPXS7v817tKlS9Svf//+Fr/66qtR24svvmjxkiVLLPbfXXQfT1QmVtoAAAAAAADkED/aAAAAAAAA5BDpUagRv1Rw5513tvj++++P2jbccMNqj+FTWyZNmmTxAw88YPHAgQOjfrvssovF6623XtTWsWNHi2fOnGlxbSwD92kDpObUjC7h7969e9S2/fbbW/zYY49ZzHL+uuHn8z777GPxiSeeaHHDhg2jfjr/evXqFbVpelTW586aFllp7wN/Xn72s59ZrOmozZo1i/otWLDA4k8//dTiQs5f0vj49FZ9b3z99ddRG8u7sys2RVDH44gjjojaTj/9dIvvu+8+i2+77baon6ZOoXo+9WyrrbayWK+br7zyStRv6tSpFq9cuTLTc9XGe2HttdeO2nRu6niTFrmqtJSYtdZay2I/nmvyeqdjzXW2cDrGrVq1itp23313i33ak37X0PtVvz1EUmpqCCEsXLjQ4sWLF1t88cUXR/0mT56c+PpRGVhpAwAAAAAAkEP8aAMAAAAAAJBD9ZoelVZ5yLfpcj+Wb+ZH48aNo8d77723xU2aNInadNy++uqrav8mhBBeeOGFav/m/fffj/r985//tHiDDTaI2rSCzaxZsyyutHSKupaWxqLS5qzOe78s9YsvvrB4+vTpFrP8t/boGLZu3TpqGzZsmMVpVd40RVBTcUKIxzdt3PQ9kvZ+4fr/I63u1a9fP4vbtm0b9bvhhhssLraCl9L3wvrrrx+16Xto3rx5UVttV/IrN1mvp2k0VXHkyJFR20YbbWRxu3btLB49enTUj/So6mkqzJlnnhm1XXjhhRbrHNNraAghPPPMM5meS98LPhVL72P0ufz7R+eir/inc/OJJ56w+PPPP8/0+sqdprPssMMOUZtWtXz++ect1vMYQvrnXVKKqb9HzXrPyudi4fS8a+XSu+++O+qn1Wj9HEu6ZqeNh763Qghh0003tVivyyNGjIj66fYTfJcpXtYxzOM5ZqUNAAAAAABADvGjDQAAAAAAQA7xow0AAAAAAEAOrfE9bbREqZZHCyGEIUOGWOzLQ2vZWM0J9vuc6N4Kuh/GihUron6aq+bzhTVvOes+O7pHi28rZ76k6z333GNx+/btozYd06FDh1qsJbnT9OjRI3qclges+91Qkjvf9H3hS37PnTvXYi2LmKbY8qiVRM+RztN777036te1a1eL9Vroz6m2XXHFFVGb7qMxadIkiz/44IOoH/O0cG3atLF4wIABFr/00ktRP92nIm0+pO2pon+nn2/NmzeP+p188skW33TTTVHbnDlzEo+P4vgx089Jfx+l81Tvj9jDJhvdp2/QoEFRm97bvvzyyxb7a2rSPglp+yz4OZs0hxs2bBg9/sMf/mCx7ocRQgjjx4+3+Omnn7bY711WKfeyIcTn/IgjjrD4tttui/pp+fSPP/7Y4v333z/qp/sppo2hvicKeR8kHa/SpJ2ztL2hdtppJ4vHjRtn8brrrpv4XP4869jpXPHfOZcsWWKxXntDCOGjjz6q9ng333xz4nNVKj/Wer3S/dy6desW9evbt6/FvXr1itp0Xz7dO9Off/0+8uWXX0Zta2rfXVbaAAAAAAAA5BA/2gAAAAAAAOTQGkmP0uVMWjrNlyDUpYVNmzaN2nbbbTeLNY1q5cqVUT/9Oy0vumzZsqiflojW1xRCvFRNUzJ8WWlNzTr77LOjtvnz51tczssW/dJZ/b9PO+20qE3Lg/vlgUm0NJ5fjqzvK3+8qVOnZjo+aiat3GgSv/xaSxW3bNkyarvmmmsszpo+kzXFo5L4ZcFbbrmlxY8//rjFmm5T3d8l0WX5nTp1itpuvPFGizVNx5fM1dTKSh2n1fGlQvUcbrbZZhaPGTMm6pd0vU1bapy1PHufPn2itn333ddiTcFAzSRd1zSdO4QQLrroIovTUl3OOOMMiwtJj8qarlGOOnToYLEvda8pUQMHDrT4008/TTyejk9aWm9aWoT+nU/P0TLG77zzTtT2pz/9yeKlS5dmOn65j7fe448aNcritHQZTf/1190jjzzSYr03DiEe09pIeyn3sfH03mSrrbaK2vRz8t1337XYf1/U73o6f/3WHfr9ceLEiVGblnmfPXt2tccOIf4Mzjqfs95Tlzu9v9x6662jtltuucXibbbZxmJNVw0hPU1O6Xupf//+UZv+HuDvbUaPHm2xjn1tz0tW2gAAAAAAAOQQP9oAAAAAAADkED/aAAAAAAAA5NAa2dNGc7p0TwPdwyCEeG+ZTTbZJGrTsnpaSm277baL+jVq1MhizRvUvNMQ4lKYfq8MLeWlecs+l04fDx48OGr77W9/a3El5Zrq/+rzR/U8Z83p1L2MNJ/cP9dTTz0VtflS5KhfOqa+LKmW1vS5pjNmzLA4bR7p3/ljUM42LlEaQgh33HGHxbqPjT93WfdV0H7ffPNN1Ka5xc2aNbP4uuuui/pNnjzZYl8OHN/bddddo8eHHnqoxVpS9KGHHor6ZS3bm3X/Ch1TvzeRftb6zwDUDh0nvw/VjjvumPh3n3zyicX//ve/Mz2XvyZUUtlZvZ8MIYTLLrssse3yyy+3WPc7TJN2fU2bfzr+rVq1svj666+P+uln7U033RS16R43lXSPqvyeT+edd57Ffs+iJDoWXbt2jdp0jl1wwQVR25///GeL066T7NFXvd69e1t86623Rm16bq+66iqL/b2J7jN04IEHWuz3Q9G/++qrr6K2pDGo5LEphv9eoPsuXnzxxRbrPl0hhNCkSZNqj6e/NYQQwqJFiyz2c6p169YW6zXB72XVpUsXiy+55JKobYsttrBY993NuhdnVqy0AQAAAAAAyCF+tAEAAAAAAMihNZIepTRV4cUXX4zaXnrpJYuzLs3WNKcQ4lJvuoxNl5CGEEL37t0tfuaZZ6I2LdelZTG7deuW+Jp8WUeWxq16DnT5b1qp6I4dO1r84IMPWuyXLGoKlF/6q+8zlpfWnaznVvu1bds26qflFZcvXx61LV68uODXlDUVpNzpOddSiCHEKRRp5Q91zmrJykcffTTqp2W9fWriyJEjLd55550tXm+99aJ+Wt79+OOPj9oqOd1RU9uGDRsWtek1UUtOppUZTpOWSqjvpx122MFin8qsy5Jfe+21ol4H0j+bdAn3scceG7U1btw48e+0PK1PFVCVVObZ0/99zz33jNr22GMPi/V6GEJcMjirYlPNNJVAl+lr+mkI8b2tpoyEUHnj+gMdX02LCCGeSzo2viS63pfodxD/maYpVn/84x+jtjfeeMPiqVOnVvu8aa89hMoaQ//evvPOOy1u2rRp1Kaff0nfO/xj/b6oW3CEUNnXw7qk4zZmzJiobf/997dYv9f7exRNC9fvgVdffXXUTz/vfEqV3g/rNf/EE0+M+mnqt0/n0s8GPb6/dtQUK20AAAAAAAByiB9tAAAAAAAAcmiNp0cpvwww61JRTX949913ozZduqax7hwdQggzZ86stl8I8VK4tKXGWmXq/vvvT3wd+F7SOfHnWJdw6/I5//7Q5W8+xY0UmbqTlBKVtdqFVosKIV5KqNWiQlh1mWqSpHlfyXRJqVY9CWHVVMMf+HOny8B1B/9777036qfpSz7d6uSTT7b42WeftdgvJT/ssMMsXrBgQdSmz11pc7tnz54Wd+rUKWp78803Lb7lllssLmSJvcpaoU0rRvkKLDp2Pt0R6bJeu3RuDxo0KGrT8fVzRSu2ZX2uSruerrXWWhYfffTRUZteN/0c06XzWu20WGn3ob/85S+rfY3Lli2L+p1//vkWU8nte3odO+aYY6I2rRajlYX22WefqJ9WYdOqQz4FTSuM+Uo0ffv2tXjatGkWZ63SWAl0rHyFqObNm1usqWYhxJWJi9kqIa1fWopa2t8l/U0l0c+tEOJ0Jq2EGUJ8LnXLktNOOy3qN2nSJIs1xS3tHOt39xDi3wP0vtQfI20rAZ3fWgn7s88+i/rVdOxZaQMAAAAAAJBD/GgDAAAAAACQQ/xoAwAAAAAAkEP1uqdNsYrZv6KQ0opaykvzhdPyw99+++3Mx0fsuOOOix5reXYdX58vfMUVV1j87bffJh6fcn31T/cJGDBgQNSmOcc+bznr/iWM66patmxpca9evRL76bl7/fXXozbdO2HWrFkWp803f62dO3euxY888ojFfr8IzXfWfXBCiEuKv//++4nPXQ58KckLL7zQYp9TrddAzef29O/0PKeNo6fvp969e1vsx1tL21ba/kM1lfWzSssMt2vXLrGf38tvzpw5mY7P9fR7a6+9dmKb3xfshhtusFj35fN7PmmJWh0fHdMQQmjTpo3FBx98cNSme1vp3H7ssceifv56jng89ttvv6hNr43jxo2z+IMPPkg8nn6m+XLEyu93ontgFLvPWDne2+r/tPXWW1v885//PPFvtOx9CHHJ77Q9Z7LuQaNzzH/e6TGKHcdK4a9xhxxyiMX+3kbvHfSeQudbCPE9jJ5/f93VfWY222yzqE2vp7pfn3+9yr8P9DcA3TOntsedlTYAAAAAAAA5xI82AAAAAAAAOVQy6VE1LTOc9ne+DNnIkSMt3njjjS32KVDDhw+3uJD0q3KWdbmmlrC8/PLLE4+hy1J//etfR/2yLkUtxyWkpUDPe+vWrS3u0KFD1E/Lo2pJ6BCSxyut7CK+p0uL/VJ+PV/Tp0+3uE+fPlE/XxrxB4Wcf12+eu6551rsl6ZryVxfHlWXRt9xxx2Jz1UOWrRoET3eYostLNbUihBCeO655yzOOgf0upn2N3558bHHHmtx06ZNLfYlLZ966qmCX1Mly/r5pOOhZYa1rHAI8b2ITyn++uuvi36dlUKvVz5dd6eddrJYS8P6tt13391iv+z/m2++sVhLdPv7UB1XX/Jb3zP6eq+99tqoX9o9UqXS+ZGW6pKUdhFCPE91XvoxTDp2CHFaY7HfH4r5XpQ3/n/Qe5URI0ZY7O8J9FqmJdNDyP7/J2214V9T1jLspXTe64OmWIcQzxd/7vTa1aNHD4vPOuusqN8222xjsc4pnwKlJeL99VTHV1+Tv3bra1y6dGnUNmrUKIu/+OKLUFdYaQMAAAAAAJBD/GgDAAAAAACQQ7lKj6qNXbizpsXosqeePXtGbSeccEK1/XT5Uwjp1ToqVdoSQ13u+/DDD1u8/vrrR/1013Bd7qtLiVcn6X1AWk3NZT1numz48MMPt9hX5HjmmWcs1l3/vaw7/eN7Xbt2tdinuujS0yOOOMLitGWdev79MnBdXppWMUjncNZqRyGE0K1bt8S+5canXSifirRy5cpMxyxm+b1/z/Tr189iHZ9XXnkl6vf5558X/FyVLOv1VNMGtKqbnyuaNuBTCZPeB3wu/kj/d00/DCFO09x3332jtiFDhli8ySabWOyvh1qF74EHHrC4SZMmUT+trKLpAZ5+fs6YMSOxn+ffNz8o91R/HY933nknattzzz0tPvvssy32aYWatj948OBMz+XnlH7v0NTghQsXJh6vHPnzoqkreg/j3696Pezbt2/UNn78eIv1/eyPUe7v9bzx721Na9thhx2iNq04e9BBB1ms96shxGOadu+p9zNZv0v494fO+4EDB0Zteu0tpCpnoVhpAwAAAAAAkEP8aAMAAAAAAJBD/GgDAAAAAACQQ7na0yZNMWXV0v5Gc9oOPfTQqE1z3z788EOL77nnnsTjY1U+b/DSSy+1WHOHfb+5c+daPGbMGIvT8hU9zUXU3Ej/XJTErDsNGza0+OCDD07sd/fdd1tcyBgXo5LKv6ftA6M5+lpyPY2eO5/rq4/9HNOcYy2L26ZNm8Tn8scfPXp0ptdYDvzeXfrY5+S3bt3a4vnz51uc9b3tx0ofd+jQIWrT0pp6fC3xHQLX1Lqiex117tw5sZ/uG6A5+Cic35vgjTfesFjvU0II4fbbb7dYP/v8Z5qW/E7bb+ORRx6x+NFHH43adO+xI4880uJC5l7a9byc6bXrb3/7W9Q2aNAgi/Xz6YYbbkg8no7nk08+GbU9/fTTFp9zzjlRm+5jc99991ncu3fvqF/aPVHa526p+vLLLzP10+9p48aNi9ouuugii9966y2L/Z5wOq46n/0+m61atbLY7yO3YMECi3U/t0qaU1n5e83+/ftb3KVLl6jtmGOOsbhXr14Wa1nvEOLPuOeff97i7t27R/30sd+TUen4/vWvf43azjzzTIv9e2lNfZ9gpQ0AAAAAAEAO8aMNAAAAAABADuUqPaqYFKi0Y6RZZ511LPbLEZcuXWrxr371K4uXL1+e+XVUUhpGknXXXTd6fMYZZ1isS4H9csM+ffpYrMsNk0pUro6e/7pOv8GPmjZtarGmWviSwM8++6zFxaQ+rk65LBteHT8/9thjj8Q2PSdaYtOX4U66jvmlv9pP0xFDiOezluL0/dS8efOix5qWUO788vCZM2da7MtiXnbZZRbrPPJLd7WcsJYj1jkaQgjvv/++xbo8OYQQmjVrZrGO/5QpU6J+xaRmVepnZBp/3dpxxx0t1vH05+6JJ56wWFM30nD+C+evgXqus5535cd7r732sthfE6655hqLfTplkrT7p0odfy2XHkIIf/nLXyw+6qijLG7UqFHiMV555ZVq/yaEEFasWGGxLy9+1113WazpHxqHsGrZ+STlMoZ6j37bbbdZPHTo0Kifprj4ewlN5U1K6w0hhMMPP9xiHStNRQ0hTsXyKYg6NzUd7rjjjov6aUpjpfLnX79Ta2pTCCG88MILFuv599ddPaaWgdd01RBC6NGjR+IxNG1LU6D+/ve/R/3y8P2RlTYAAAAAAAA5xI82AAAAAAAAOcSPNgAAAAAAADmUqz1tstLcX58HnLQvjs951FzJTp06RW0TJ060+KGHHqr2eGmvKYQ4f1jz58ol7zSLgQMHRo917wylecQhxGXWa+N8UYK2ZtLmmPI585qT37x5c4tffPHFqN+SJUtq+AoRQpz3G0IILVq0yPR3fk8plXX+6XP7eX/jjTdarHtxeJpTPmDAgKhNS5SXO79/hX5WnXrqqVFb165dLdZy6prb7R/rPNX920II4b333rN48803j9r0OqC53XPmzKnmv6ieHkP3JPDX6HL+nMx6PfX3LBdffLHFOt/8/inXXnvtao9dE+xFVJi0PdV0HAcPHhy1XXXVVRb7vcYmT55scVpp4bTnpiTxqudVr6963e3Xr1/iMa677jqL0/a91O8SIcTXWt1nTMtVhxDCIYccYnExeyWVGv1suf766y1u2bJl1O8Xv/iFxWuvvXbUpu/7tL2cdG/TtHsT5e+zdL8jLWGt++WEEO9hxNxbPT1HaedLx/f444+32N+Haj8/j4YNG2ax7mOThz1sPFbaAAAAAAAA5BA/2gAAAAAAAORQgwJL5+ZiLWzakk/9f3Q5lF8qNXbsWIv90uy+fftanLUccdYlz8WqqqqqlbrFdT2GWub7rbfeito22mgji3W5W8+ePaN+Pn3mB2nLHOtiuaGOqf5fvkzuwoULE1+HS5ObUVVV1SPUgroex6xzTPnl/E8++aTFOsaXXHJJ1O/qq6+2OK2UdJbXUJ3aXs6f17nol/dq2lnDhg2jNi0/2aZNG4v98m59/+oxdC6HEMKIESMsPvjgg6O2pHKpfqzvu+8+i325zDpIcSyZuajnXdMMQ4iX1W+55ZYW+2W9Oq5a3nLRokVRPx2rCRMmRG2aiqUpXO3bt4/6ffLJJ9X8F6vS/8u/3qzLkvM6F6s5fmJb0jVpgw02iB7Pnz/fYp3rM2bMiPrttttuFqelU+QotSm3c9GPW9bPkqz9dA7Mnj07attiiy0snj59etSmY/ztt99meh1eud2j6meV/t9pJYKree7Vxv4YSXHa6wshhP3228/iMWPGJP6d9tPy4qt7viLldi76exi9D/f3B3ovqqWer7zyyqjftttua7GOT1ppav86ku5vpk2bFj3eddddLa7rtJv6not1yc9FTZO74447LPbjpF577bXo8S677GLxZ599VtOXWFuqnYustAEAAAAAAMghfrQBAAAAAADIoZKsHpV1SaDuNn7TTTdFbboU7tVXX43aND0n63PlaKnxGuWXqukys/XXXz/x73R54I477hi16RJQ7bfeeuslHs+ndSSlU/hd5tu1a2fxNttsE7UdcMABFmvKnK/Mcvnll1s8adKkqK2QyiqlTlPIQghhu+22s1jfJ//5z3+ifnU9dyplbvr/Uytj+KWiunz4wAMPtHjmzJlRv1NOOcVirei04YYbJh7PV1dIeo2ffvpp1DZ8+HCL87hrf33RFJfFixdHbfpYx67Yzy2dp/4z8+abb672GGnLkNPoNTot9bXcZB0b/3mkn116DJ9aoe+XtPGlCtTqpaVHpcmaMqMV1Hw6nI7jSSedFLUlpUSlvd5yqFijny3+MyLpc6fY/1vHLWuFLh3PEOLX6I/x1FNPWfyPf/zD4qOPPjrqpyk9vs1XvypnPtUzaxru448/brFPu9b0KB1vTR8PIYQpU6ZYvPvuuyc+lx7jo48+SmzDqrJuMeLvPbV6W9q9iN4rHXnkkVHbsmXLMr/O+lY5d0oAAAAAAAAlhB9tAAAAAAAAcogfbQAAAAAAAHKoZPa0yZp/rXsr3HrrrRanlcM96qijorZKyhOtbVoK1uegar6vxqNGjYr6/e53v7N4xYoVFvtS2/qe0DK2IcT74ujfdenSJeqn++74ktVZ89eHDh1qcceOHaO2IUOGZDpGnmXNxfW5vjrnNAd/7ty5mY+vbVnHo1L569a4ceMsPv3006M2zf0dO3asxX5vEd0noNjzr2OoucMXXnhh1E/LMJL/Xb2sc6U2ju/3zdDx12tlsfvR6B4P5bDfRpqsY6Pn2Ofd63nW/YAmTpxY1HNh9Yqdb1nLSnfu3Nlif//x4YcfWjxv3rzU11nM6yhFukdM2ueRzo9CzkHSMbOOdSHPpffHDz74oMWHHHJI1K9bt24Wb7311lFbMfuYVRp9zzz99NNRm+5b2bhxY4v9XNRy3X5fzKTn8nvClftnXE1lvWaeccYZUZvudarH8PucHnbYYRb7kt+lNHdYaQMAAAAAAJBD/GgDAAAAAACQQyWTHpXEL80eNGiQxXvvvbfFmmYTQrwE0S89LaWlUvXNn6tZs2ZZrGUMQwjh0EMPtVjLZvsyba1atar2+GnLYX158c0339xiXSrrlz2mlZDUx7qU9euvv476aQlrLZFcyrLOAT1/p556amI/Le+8aNGiOn1Nlcovv73oooss7tWrV9TWvXt3i4st2Zz03CtXrozadH6cfPLJFmsKYwjxPEXNpJX+TVvOr/3089O36Rj7z1YUT9OGe/bsmdhv6dKlFvvSsmm4hhbGn69izp+fi3rvo5+ZaddhTQEIIU6F1dfk74eLTd0pBf7/SSqDXoiapmD7e8g0+vqfe+45iydNmhT10/Qon4L+6quvWlxp2zroez1tnmqsKeMhxOlmJ554osU6R6t7rPTe5+GHH7b42WefTfwbFEa/E5511llRW1La8AUXXBD1mz59usWFzNO8YaUNAAAAAABADvGjDQAAAAAAQA7xow0AAAAAAEAOlcyeNkl7m7Rv3z7qd+2111qsJdw01zCEEKZNm1btsVEzX3zxhcWDBw+O2v71r39ZfM4551jsS+gtXrzY4o022sjidu3aRf10PxW/j4aWdNMc4XXXXTfqp+Wnn3zyyahNcyVbt25t8dSpU6N+bdu2tbjS3kt6PjfddNOoTcdk/PjxFvs9gVA3vvzyS4v79OkTtT3wwAMWa5687qmRxu8/o/tq+FLe999/v8W6/0mlzZWsaqO0fdox0s67jv9mm20WtWnu/ttvv20xe9rUHt1zrXnz5lGbnv+FCxdazF5QpaVZs2YW77vvvhb7/fb0nsPvzzBixAiLdb84v1dDbezzUknq6zNJ75t///vfR23Dhw+32H+OT5482eLZs2dbXMp7diTxn2lNmza1WPecDCF5zyf/WTV06FCL+/fvb/HGG2+c+DqWLFkSPZ44caLFut8K97k1o9fDBx980GIdd++tt96y+K677ora/HukVLHSBgAAAAAAIIf40QYAAAAAACCHSiY9SmA7+LYAAARySURBVOkyuUsvvTRq0yXFupx47NixUb9yXD6YB2lLEe+8885q46x8OcsmTZpY7JcB13b5Qy1b7L3++uu1+lylRFPUfKlKHZ9hw4ZZzNxb85YvXx493n///S3eaaedLD7ppJOifltttZXFmu6m6YchhDBy5EiLdYlqCKuWIke6tDLcWZfvF7vMX6+xH374YdSmy8Lvvfdei0nBqD1aWnbZsmVRmy7vfu+99yz2aTWMR774uajX4kWLFlmsadYhxHOxd+/eUZumcmvaOalyNZO0DYNPzantNCo93vz586M2LQfutxwYMmSIxeeff77F/vO+HPhzrulH/h4j6/h89tlnFvfr18/igw46KOo3ZcoUi2fOnBm16XcN7m1rT4cOHSzebrvtEvvp+0Cvk7o9QDlhpQ0AAAAAAEAO8aMNAAAAAABADjUoZJlfgwYNclHuQ3f2fuONN6I2rUSky/k7deoU9VuwYEGm5ypmaXpdqKqqqnlJkZCfMaxQM6qqqnrUxoHyMo66hLtRo0ZRmy7TL6cl+8zFul8uvgaU3VysDZruuP3220dtPXv2tHjChAkW+zSqNfleKLe5qNW7tLJQCCEccMABFo8ePdpirRoTQkku0a+ouajXzs6dO1t85plnRv30s1XTEUMI4eWXX7ZYUzy4Ry34uaLHSdUTs96/+BT+2kgNbtGihcWjRo2K2ubNm2fxY489ZrGm8xSoouZiuSrFudiwYcPo8bvvvmtxy5YtE//u9ttvt/iUU06xuATvSb1q5yIrbQAAAAAAAHKIH20AAAAAAAByiB9tAAAAAAAAcqhk9rTRXPvf/OY3Fl955ZVRP81RnTNnjsXdu3eP+qXlqCaV+qvP0rWlmKOIVZAvXAaYi2WBubgafr8HlZd88XKei6Vw/mtJxc7FrGWlS2G88zoX0+ZRmmLOud/TJu0Yun9O1u8jHTt2jNr22msvi8eOHWuxlkEuUMXOxXKS17mYpnHjxtFj3a9W97FdtmxZ1K9t27YWr1ixoo5eXb1gTxsAAAAAAIBSwY82AAAAAAAAOVR9fbsc0mWGWubOl7fUZYbnnnuuxd99911Rz1UKy1IBAKhNfPbVL85/+eNes+6tyfNayBYKWcuI6+t/8803o7aPP/7Y4hqkRAH17ptvvokejx8/3uJjjjnGYv1eH0LZpUStFittAAAAAAAAcogfbQAAAAAAAHKIH20AAAAAAAByqGRKfifRsnleIfvY5F0plnDDKiinWAaYi2WBuVgGmItlgblYBpiLZYG5WAaYi2WBkt8AAAAAAAClgh9tAAAAAAAAcqjQkt+fhBDeqYsXUqxySoFKsWktHit3Y1hBGMfSxxiWB8ax9DGG5YFxLH2MYXlgHEsfY1geqh3Hgva0AQAAAAAAwJpBehQAAAAAAEAO8aMNAAAAAABADvGjDQAAAAAAQA7xow0AAAAAAEAO8aMNAAAAAABADvGjDQAAAAAAQA7xow0AAAAAAEAO8aMNAAAAAABADvGjDQAAAAAAQA79PwW4VVyfmyBUAAAAAElFTkSuQmCC\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: 1min 5s\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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0018, PSNR (dB) = 27.44\n",
      "MSE = 0.0045, PSNR (dB) = 23.46\n",
      "MSE = 0.0010, PSNR (dB) = 30.11\n",
      "MSE = 0.0020, PSNR (dB) = 27.01\n",
      "MSE = 0.0035, PSNR (dB) = 24.51\n",
      "MSE = 0.0006, PSNR (dB) = 32.38\n",
      "MSE = 0.0050, PSNR (dB) = 22.97\n",
      "MSE = 0.0032, PSNR (dB) = 24.95\n",
      "MSE = 0.0049, PSNR (dB) = 23.06\n",
      "MSE = 0.0033, PSNR (dB) = 24.77\n",
      "Wall time: 18 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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9643\n",
      "0.9590\n",
      "0.9641\n",
      "0.9813\n",
      "0.9398\n",
      "0.9793\n",
      "0.9160\n",
      "0.9430\n",
      "0.9432\n",
      "0.9698\n",
      "Wall time: 3.39 s\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": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9970\n",
      "0.9942\n",
      "0.9981\n",
      "0.9972\n",
      "0.9932\n",
      "0.9990\n",
      "0.9906\n",
      "0.9943\n",
      "0.9876\n",
      "0.9913\n",
      "Wall time: 13.1 s\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": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Simple Autoencoders in Keras - 784-64-784-tf1-rpi.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
}