Deep Autoencoders in Keras - 784-256-128-64-128-256-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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using plaidml.keras.backend backend.\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "os.environ['KERAS_BACKEND'] = 'plaidml.keras.backend' # https://github.com/plaidml/plaidml/blob/master/docs/index.md\n",
    "\n",
    "from keras.models import Model\n",
    "from keras.layers import Input, Dense\n",
    "from 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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:plaidml:Opening device \"opencl_intel_hd_graphics_4000.0\"\n"
     ]
    }
   ],
   "source": [
    "input_img = Input(shape=(784,))\n",
    "encoded = Dense(256, activation='relu')(input_img) # apply to previous layer = input_img\n",
    "encoded = Dense(128, activation='relu')(encoded)\n",
    "encoded = Dense(64, activation='relu')(encoded) # the bottleneck = encoding_dim\n",
    "\n",
    "decoded = Dense(128, activation='relu')(encoded)\n",
    "decoded = Dense(256, activation='relu')(decoded)\n",
    "decoded = Dense(784, activation='sigmoid')(decoded)\n",
    "\n",
    "autoencoder = Model(input_img, decoded)\n",
    "autoencoder.run_eagerly = True # plaidml / keras not running eager mode by default\n",
    "\n",
    "encoder = Model(input_img, encoded)\n",
    "\n",
    "encoded_input = Input(shape=(encoding_dim,))\n",
    "\n",
    "decoder_layer = autoencoder.layers[-3](encoded_input)\n",
    "decoder_layer = autoencoder.layers[-2](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-1](decoder_layer) # extract last layer\n",
    "\n",
    "decoder = Model(encoded_input, decoder_layer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 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": [
      "1595253818\n",
      "Train on 60000 samples, validate on 10000 samples\n",
      "Epoch 1/50\n",
      "60000/60000 [==============================] - 26s 432us/step - loss: 0.2142 - val_loss: 0.1389\n",
      "Epoch 2/50\n",
      "60000/60000 [==============================] - 21s 356us/step - loss: 0.1236 - val_loss: 0.1109\n",
      "Epoch 3/50\n",
      "60000/60000 [==============================] - 20s 338us/step - loss: 0.1070 - val_loss: 0.1018\n",
      "Epoch 4/50\n",
      "60000/60000 [==============================] - 20s 336us/step - loss: 0.0999 - val_loss: 0.0964\n",
      "Epoch 5/50\n",
      "60000/60000 [==============================] - 20s 338us/step - loss: 0.0953 - val_loss: 0.0921\n",
      "Epoch 6/50\n",
      "60000/60000 [==============================] - 18s 295us/step - loss: 0.0920 - val_loss: 0.0900\n",
      "Epoch 7/50\n",
      "60000/60000 [==============================] - 20s 338us/step - loss: 0.0897 - val_loss: 0.0878\n",
      "Epoch 8/50\n",
      "60000/60000 [==============================] - 20s 339us/step - loss: 0.0879 - val_loss: 0.0865\n",
      "Epoch 9/50\n",
      "60000/60000 [==============================] - 20s 338us/step - loss: 0.0864 - val_loss: 0.0850\n",
      "Epoch 10/50\n",
      "60000/60000 [==============================] - 20s 340us/step - loss: 0.0853 - val_loss: 0.0839\n",
      "Epoch 11/50\n",
      "60000/60000 [==============================] - 20s 339us/step - loss: 0.0840 - val_loss: 0.0829\n",
      "Epoch 12/50\n",
      "60000/60000 [==============================] - 21s 342us/step - loss: 0.0832 - val_loss: 0.0821\n",
      "Epoch 13/50\n",
      "60000/60000 [==============================] - 20s 341us/step - loss: 0.0823 - val_loss: 0.0812\n",
      "Epoch 14/50\n",
      "60000/60000 [==============================] - 21s 353us/step - loss: 0.0815 - val_loss: 0.0764\n",
      "Epoch 15/50\n",
      "60000/60000 [==============================] - 22s 368us/step - loss: 0.0807 - val_loss: 0.0804\n",
      "Epoch 16/50\n",
      "60000/60000 [==============================] - 21s 348us/step - loss: 0.0797 - val_loss: 0.0791\n",
      "Epoch 17/50\n",
      "60000/60000 [==============================] - 22s 361us/step - loss: 0.0794 - val_loss: 0.0785\n",
      "Epoch 18/50\n",
      "60000/60000 [==============================] - 21s 356us/step - loss: 0.0788 - val_loss: 0.0781\n",
      "Epoch 19/50\n",
      "60000/60000 [==============================] - 22s 366us/step - loss: 0.0784 - val_loss: 0.0777\n",
      "Epoch 20/50\n",
      "60000/60000 [==============================] - 21s 355us/step - loss: 0.0779 - val_loss: 0.0775\n",
      "Epoch 21/50\n",
      "60000/60000 [==============================] - 21s 355us/step - loss: 0.0775 - val_loss: 0.0770\n",
      "Epoch 22/50\n",
      "60000/60000 [==============================] - 22s 371us/step - loss: 0.0771 - val_loss: 0.0766\n",
      "Epoch 23/50\n",
      "60000/60000 [==============================] - 21s 352us/step - loss: 0.0767 - val_loss: 0.0768\n",
      "Epoch 24/50\n",
      "60000/60000 [==============================] - 21s 352us/step - loss: 0.0765 - val_loss: 0.0759\n",
      "Epoch 25/50\n",
      "60000/60000 [==============================] - 20s 330us/step - loss: 0.0761 - val_loss: 0.0758\n",
      "Epoch 26/50\n",
      "60000/60000 [==============================] - 21s 347us/step - loss: 0.0759 - val_loss: 0.0760\n",
      "Epoch 27/50\n",
      "60000/60000 [==============================] - 22s 359us/step - loss: 0.0756 - val_loss: 0.0755\n",
      "Epoch 28/50\n",
      "60000/60000 [==============================] - 19s 316us/step - loss: 0.0754 - val_loss: 0.0748\n",
      "Epoch 29/50\n",
      "60000/60000 [==============================] - 20s 330us/step - loss: 0.0751 - val_loss: 0.0749\n",
      "Epoch 30/50\n",
      "60000/60000 [==============================] - 22s 371us/step - loss: 0.0746 - val_loss: 0.0747\n",
      "Epoch 31/50\n",
      "60000/60000 [==============================] - 19s 315us/step - loss: 0.0747 - val_loss: 0.0744\n",
      "Epoch 32/50\n",
      "60000/60000 [==============================] - 20s 336us/step - loss: 0.0745 - val_loss: 0.0746\n",
      "Epoch 33/50\n",
      "60000/60000 [==============================] - 20s 341us/step - loss: 0.0744 - val_loss: 0.0744\n",
      "Epoch 34/50\n",
      "60000/60000 [==============================] - 20s 330us/step - loss: 0.0742 - val_loss: 0.0737\n",
      "Epoch 35/50\n",
      "60000/60000 [==============================] - 20s 335us/step - loss: 0.0740 - val_loss: 0.0740\n",
      "Epoch 36/50\n",
      "60000/60000 [==============================] - 20s 337us/step - loss: 0.0738 - val_loss: 0.0738\n",
      "Epoch 37/50\n",
      "60000/60000 [==============================] - 20s 339us/step - loss: 0.0737 - val_loss: 0.0734\n",
      "Epoch 38/50\n",
      "60000/60000 [==============================] - 20s 336us/step - loss: 0.0735 - val_loss: 0.0734\n",
      "Epoch 39/50\n",
      "60000/60000 [==============================] - 20s 336us/step - loss: 0.0734 - val_loss: 0.0735\n",
      "Epoch 40/50\n",
      "60000/60000 [==============================] - 20s 338us/step - loss: 0.0733 - val_loss: 0.0731\n",
      "Epoch 41/50\n",
      "60000/60000 [==============================] - 20s 340us/step - loss: 0.0732 - val_loss: 0.0729\n",
      "Epoch 42/50\n",
      "60000/60000 [==============================] - 20s 336us/step - loss: 0.0730 - val_loss: 0.0730\n",
      "Epoch 43/50\n",
      "60000/60000 [==============================] - 20s 337us/step - loss: 0.0729 - val_loss: 0.0727\n",
      "Epoch 44/50\n",
      "60000/60000 [==============================] - 20s 337us/step - loss: 0.0728 - val_loss: 0.0727\n",
      "Epoch 45/50\n",
      "60000/60000 [==============================] - 20s 337us/step - loss: 0.0727 - val_loss: 0.0728\n",
      "Epoch 46/50\n",
      "60000/60000 [==============================] - 20s 338us/step - loss: 0.0725 - val_loss: 0.0724\n",
      "Epoch 47/50\n",
      "60000/60000 [==============================] - 22s 373us/step - loss: 0.0725 - val_loss: 0.0724\n",
      "Epoch 48/50\n",
      "60000/60000 [==============================] - 21s 355us/step - loss: 0.0724 - val_loss: 0.0709\n",
      "Epoch 49/50\n",
      "60000/60000 [==============================] - 22s 368us/step - loss: 0.0723 - val_loss: 0.0722\n",
      "Epoch 50/50\n",
      "60000/60000 [==============================] - 21s 349us/step - loss: 0.0722 - val_loss: 0.0724\n",
      "Wall time: 17min 23s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<keras.callbacks.History at 0x21016a98f88>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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+WH4yJAAATX0lEQVR4nO3b+XON5//H8evkJJFII5FFQtHREqU6jCqGGtWqWgZDp7WvJbSKWIpQ26C0I6UMxpqIpmrEULRqSkNrq9rbalFTFbGklizE7nz/gK/kc71mTpqrnefjx87TO0fuXOe+z1vq8fl8BgAAAAAAAG4JKOsXAAAAAAAAgP+PpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgIMClTg8PNwXExNj3RcWFkov5s6dO1IfEhJi3V67dk2aXaNGDamPiIiQ+tzcXOs2Ly/P3Lp1yyN9gWIEBAT4vF6vdR8fH6/Ol/obN25Yt+rPU2hoqNTfvn1b6itVqiT1ubm5V30+X6z0h4oRHh7ui46Otu6Vc2uMMSdOnJD6ChUqWLfqWVS/zx6PdlSUPj8/3xQVFfnlLIaFhfmioqKse5/PJ81Xv2/nzp2zbpVza4wxcXFxUh8bqx2T/Px8qc/OzvbbWQwJCfGFhYVZ90FBQdJ85WwZY0xRUZF1m5OTI81OSEiQevWerlz3c+fOmatXr/rlLJb2++lvv/0m9ZGRkdbtxYsXpdnVq1eXetWjR4+k/sKFC347i4GBgb7g4GDrXjm3xujvwcr7+5kzZ6TZtWvXlvp79+6VWn/jxg2/PaOGhob6lPc89ZlTvS+eOnXKur17926pvpZy5cpJvfqsVVRU5LezqD7fqO8bFStWlPq8vDzrVr0vqp+T1M8myvfRn/fFsLAwn/J9Vq+hel6U91/1GVX9zK9cE2OMOXv2rNTn5eU99ixKS5uYmBgzbdo063737t3KePnBRrlprV69Wpr94YcfSv3rr78u9YsWLbJulyxZIs0uidfrlR44x4wZI81Xbyrr16+3brOysqTZ6geM48ePS33Pnj2lfv78+X9Jf6AE0dHRZtKkSdb94MGDpfnqg33r1q2t29TUVGl2r169pF5ZShqj/cyuXLlSml2SqKgoM3LkSOv+/v370vzhw4dL/aBBg6zbdevWSbP79esn9UOGDJH6rVu3Sv3IkSP9dhbDwsJM+/btrXt1gdWmTRupP3r0qHU7YcIEabZ6L1I/iCrXvVGjRtLskkRHR5sPPvjAulfOijHGNGnSROo7dOhg3U6dOlWaPXHiRKl/+PCh1KsP4qNHj/bbWQwODpbu+02bNpXmq4uP7t27W7fqM+SKFSukXlnKG6MtAxcsWCDNLkmFChWke76ypDPGmFGjRkl98+bNrVv1g5ny82GMMTVr1pT69PR0qT906JDfzmJUVJRJSkqy7tV/NO3atavUb9myxbpV74vq8029evWkvnfv3tatP++LFStWNCNGjLDulX8wMsaY06dPS73yDJyZmSnNnjlzptSrn/+6dOki9Zs2bXrsWeR/jwIAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwkMfn89nHHo99bIxZuHCh9GJiY2Olfvfu3dZtz549pdktWrSQ+pEjR0p9vXr1rNuZM2eac+fOeaQvUAz1Gp47d06af/DgQakPDg62brdv3y7NXrJkidSPGjVK6gMCtJ1nSkrKYZ/P10j6Q8VQr+OCBQuk+XFxcVIfEhJi3Z48eVKanZycLPUdOnSQ+sGDB1u3Y8aMMX/88UeZnMXs7GxpflpamtTfvXvXun3uueek2T169JD6lStXSv3FixelfvLkyWV2Fps1aybNb9mypdSvWrXKut2yZYs0u3HjxlI/Z84cqc/Ly7Nu09LSzKVLl8rkLK5du1aa/+DBA6mPiYkptdkdO3aU+jfeeEPq1fffgQMHltlZnDt3rjS/WrVqUp+Tk2Pd3rt3T5o9YcIEqZ8xY4bUFxQUWLfp6enm8uXLZXIWlc8wxhjz5ptvSn2tWrWs2wYNGkizu3XrJvUDBgyQ+vj4eKmfPXt2mZ3FmjVrSvPbt28v9crzutfrlWanpKRIfUJCgtRXrFjRuv3ll1/MzZs3y+QsdurUSZo/bdo0qf/222+tW/U5S/3M/+mnn0r9sWPHpD41NfWxZ5HftAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcFCgEnu9XhMREWHdr1u3Tnoxe/bskfr09HTrtkWLFtLsSZMmSf2sWbOkvkePHtZtYWGhNNufkpOTpX7t2rVSv2XLFut2yZIl0uxGjRpJ/bx586Te5/NJfUpKitT7k8fjkfpu3bpJ/fLly61b9WdqwIABUp+amir1SUlJ1m1AgP/23B6PxwQFBVn3ffv2leZnZWVJ/cmTJ63bunXrSrOvXr0q9TExMVI/e/ZsqfenoKAgExcXZ923a9dOmj958uRS6xs3bizNVt6vjTGmY8eOUt+vXz/r9v79+9JsfypfvrzUd+7cWeoXLFhg3Y4YMUKa3aVLF6nfsGGD1L/77rtSX5auXLki9WPHjpV65V6nvod9/PHHUj9u3DipV+67Dx48kGb70/Tp06U+MzNT6rt3727dqtfwo48+kvrx48dL/enTp6Xen/dRr9drnnjiCet+6tSp0vw+ffpIvfKZTv0816pVK6lXn83Wr19v3ao/IyVRn1HVz6oNGzaU+mHDhlm36vdh2bJlUp+YmCj1LVu2lPri8Js2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOClTi8uXLm4YNG1r3R44ckV7Mzp07pX7Tpk3W7ZQpU6TZdevWlXpV27Ztrdv9+/f77etWrlzZJCYmWvc+n0+ar1wTY4z5+++/rds5c+ZIs/Py8qT+0KFDUt+7d2+p96ewsDDz/PPPW/cXLlyQ5mdmZkr93r17rdtVq1ZJs69fvy71qq+//tq6zc/P99vXjYyMNK1bt7bun3zySWl+gwYNpH7GjBml0hpjzLZt26ReFRoaWqrzSxIeHm5atWpl3bdv316a//LLL0v9jh07rNuff/5Zmq2+v6u+++4767agoMBvXzc0NNQkJCRY9+qzyjfffCP127dvt26XLVsmzS5XrpzUb9y4UeoPHjwo9f5UpUoVM2TIEOs+NjZWmq9+L1555RXrtlmzZtLsWbNmSb3qqaeesm6Dg4P99nVr1Kgh/d1++OEHaf7KlSulfs+ePdbt5cuXpdnHjh2TelWnTp1KdX5JIiIipHvd4cOHpfnKM6cxxuzbt8+6XbNmjTRblZWVJfXK3/XmzZvqyylWbGys6dGjh3Wv3hffe+89qa9Tp451q17DwsJCqVdNnDhR6nfv3v3Y/85v2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAgwKV+OHDh6agoMC6v3btmvRizp8/L/VhYWHWrcfjkWZ3795d6ps2bSr1J0+etG7v3LkjzS7JpUuXzPTp0637efPmSfPnz58v9a+++qp1+9NPP0mzN2/eLPUjRoyQ+qioKKn3p/DwcNOyZUvrvm3bttL8rVu3Sv3Zs2et29zcXGl2RkaG1K9YsULqr169at0GBQVJs0tSrlw58/TTT1v3x48fl+YPHTpU6jt37mzdvvLKK9LsrKwsqU9KSpL61NRUqfenwMBAExkZad1fv35dml+rVi2pDwkJsW5nzZolzf7iiy+kfu7cuVK/c+dO6zY/P1+aXZLY2FiTmJho3bdq1Uqav3r1aqk/c+aMdZuXlyfNVs/KsmXLpL5atWpS708PHjyQvh9TpkyR5h89elTqe/XqZd0q7yHGGHPgwAGpV39my5UrZ90GBPjv339zc3PNokWLrHv1PeaZZ56ReuXsnjp1SpqtPpetWbNG6oODg6W+W7duUl+SkJAQU7t2beve6/VK8/fs2SP1yn10w4YN0mz1GVW51xhjzM2bN63bR48eSbNL4vV6pfelyZMnS/ObNGki9du2bbNuJ06cKM3Ozs6W+hkzZki9ep8uDr9pAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgoEAljomJMQMHDrTu27VrJ72YXr16Sf2NGzes2/z8fGm2ql69elIfGRlp3Xq9XvXlFKty5cpm0KBB1n1goPQjYtq2bSv1P/74o3Xbp08fafbmzZulvlq1alL/1VdfSb0/eb1eExERYd0nJSVJ88eNGyf1+/fvt25fe+01aXZGRobUp6WlSX18fLx1W1RUJM0uicfjMaGhodb9sGHDpPn79u2T+pEjR1q3n332mTS7RYsWUq++5ylnwd88Ho8JDg627hcvXizNf/bZZ6X+9u3b1u2VK1ek2apff/1V6ps2bWrdHj9+XH05xQoICDDly5e37hMTE6X5Q4YMkfodO3ZYt6NHj5Zmp6amSv3nn38u9TVq1JB6f7p79645ffq0dd+8eXNpfkJCgtQHBQVZtzk5OdJs1UsvvST1yrPZrVu31JdTrAcPHkjvS+pzWNWqVaVeuddNnDhRmq2qVKmS1K9ataqUXon/nThxQuqVz1DGGOlZq3HjxtJs9Rm1sLBQ6v/66y/rVrn//y+3bt2S3geU77Exxhw4cEDq4+LirNsNGzZIs9Vr/uWXX0p9/fr1pb44/KYNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDPD6fzz72eOxjY8zQoUOlF/PCCy9I/Y4dO6zbM2fOSLOPHDki9aXN5/N5/DFHvYa7du2S5i9atEjqQ0JCrNupU6dKs2vWrCn1/4DDPp+vkT8Gqddx+PDh0vzc3FypT05Otm4XL14szV62bJnUl7ayOotdunTxx5ct1saNG63bESNGSLMXLFigvpzSVmZncc6cOdJ89b6YkZFh3UZHR0uzU1JSpL60ldVZfOutt6T5jRppP2r9+/e3bmfPni3NnjdvntT/A8rsLCYlJUnz33//fakfO3asdfviiy9Ks0ePHi31pa2szmK7du2k+W3atJF65T1S+UxijDHp6elS/w8os7O4dOlSaX6dOnWkXjm76jNn/fr1pb60ldVZrF69ujR/+vTpUn/s2DHrtqioSJq9fPlyqf8HPPYs8ps2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOClTiqlWrmlGjRln3KSkp0osJCwuT+rZt21q3u3btkmZ37NhR6qtXry71v//+u3V78OBBaXZJqlSpYt555x3r/sCBA9L8vLw8qa9du7Z12717d2n2woULpT4jI0Pq+/XrJ/XK9/1/qVSpkunRo4d1f/jwYWn+qlWrpP78+fPW7cmTJ6XZy5cvl/odO3ZIfZ06dazbpUuXSrNLEh8fb95++23r/uLFi9L8goICqZ83b551u3//fmn2gAEDpP7OnTtS36RJE6lPSkqS+pJER0dL94usrCxp/tmzZ6U+PDzcut27d680+/vvv5f6zMxMqR82bJh127VrV2l2SeLi4kzfvn2te/UsDhw4UOpzcnKs2yNHjkizBw8eLPWRkZFSHxUVJfXJyclSX5L4+HjTv39/6149W5988onU375927o9c+aMNDs9PV3qs7Ozpf7ChQvWrXrOS6I+owYHB0vzq1WrJvUVK1a0bn0+nzR7w4YNUq9+n9X7bps2baS+JEFBQSY2Nta6V3+e1b9bp06drNu0tDRptvrMefr0aalXvo/jx4+XZpekcuXKJjEx0br/888/pfn379+XeuUe3adPH2m2+vkyISFB6gsLC6W+uOdxftMGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBHp/PZx97PH8bY/4qvZeDYjzl8/li/TGIa1imuI7/flzD/wau478f1/C/gev478c1/G/gOv77cQ3/Gx57HaWlDQAAAAAAAP4Z/O9RAAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA76Px76IimdPlCOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd7hU1dXH8UUSRbBgKGpUigjSbCgqhGY3SlQkogJGgl2JTxTFB2OLGgJWsERsGONjN0iiiSKiQnyMMYodEbtgQUBFELGg9/3D1+VvL+4Ml8vMvWdmvp+/1rg3c4+zZ59z5jx77dWgqqrKAAAAAAAAkC0/qu8DAAAAAAAAwMp4aAMAAAAAAJBBPLQBAAAAAADIIB7aAAAAAAAAZBAPbQAAAAAAADKIhzYAAAAAAAAZ9JPV6dygQQPqg9eTqqqqBoV4H8awXi2qqqpqUYg3YhzrD3OxLDAXywBzsSwwF8sAc7EsMBfLAHOxLFQ7F1lpA9Sdd+r7AACYGXMRyArmIpANzEUgG6qdizy0AQAAAAAAyCAe2gAAAAAAAGQQD20AAAAAAAAyiIc2AAAAAAAAGcRDGwAAAAAAgAzioQ0AAAAAAEAG8dAGAAAAAAAgg35S3weAyrHOOut43KtXL4932GGHpN9RRx3lcZMmTTxu0KBB0u/TTz/1eNy4cUnbxIkTPf7qq69qecSob3HMv1dVVVXHRwKlc9nM7Mc//rHHy5Ytq+vDAQAAyKQf/eiHNRLffvttPR4JShkrbQAAAAAAADKIhzYAAAAAAAAZxEMbAAAAAACADGJPGxRNq1atkteTJk3yeKuttvK4cePGOd9D80Dj/ibNmzf3+MQTT0zaJk+e7PH8+fNreMSob+3atUte33LLLR5ffvnlHt9xxx1JP/a4KY7111/f45EjR3p8xBFHJP10f6k999wzaVu4cGGRjg65xHMl8wNm6d5Tuq8C34/C0vnHZwtUBv290rZt26RtyZIl1cZm6b6b+h7x3PHNN98U5DhRulhpAwAAAAAAkEE8tAEAAAAAAMgg0qNQUFqiW1OUzMy22247j3X58Jdffpn0e/bZZz2eMWOGx9tss03Sr3v37h5vvPHGSVvXrl09fuCBB2p07Kgf+l349a9/nbS1aNHCY5bzF99aa62VvB41apTHv/3tbz1eb731kn4bbrihxzvuuGPSNmXKlEIeIoSOV6dOnTyOKaezZ8/2WJdmF2IexfLv+jouA6fUafHp+XTo0KFJ26mnnuqxppiOGTMm6cc4rZqmmpmZDRkypNr49ttvT/rp63jvU0xxnjZr1sxjnadLly6ts2MqVTH99Hvcl5SHRo0aJa/1nubYY49N2nbaaSeP9X61YcOGSb+vv/7a47feeitp0zR0/XfxN9SIESM85hxdmVhpAwAAAAAAkEE8tAEAAAAAAMigOk+P0mWFcSm+LjfVpWRm6a7ZLEHMrv32289jrRAVLV682ONDDz00aXvsscc81rFu37590u+SSy7xuFu3bknbPvvs47GmZ/DdyR49D+y+++5J22effeaxpsqhcLRawW677Za0HX300R6vu+66Od9Dx7Bnz55J28MPP+xxPK9jzeiYDB8+3GOtRmFmduGFF3ocU5ZqQ78z8dyr5+m77747adP5jOLQ74ReI83MfvrTn3qsczumR2HVBg4cmLweN26cxz/5yQ+31prubVa3aQ2aNq6pcWZmW2+9tcd/+tOfPH788ceTfqV4z1SI6l1rr722x9tuu23S1qVLF4/1vmTu3LlJP1JYsiWmtWka1OGHH+6xXi/N0m0fcqXGmaXjHcdef9/q3IttK1as8Djej2laeiGu4yg9rLQBAAAAAADIIB7aAAAAAAAAZBAPbQAAAAAAADKozve02WCDDTyeMGFC0taxY0eP494HH3/8scePPvqox9OmTUv6vfbaax5rOcV8uaW1zXnV3EbNQ6wkMb9T8yxfeeWVpE33vTjqqKM8zpfzrbmeWhbPzGyzzTbzOJaz1Lz+UszJriSaL6yxmdlzzz3n8YIFC+rsmCpJmzZtPL7iiiuSNi3lnW8e6R4OJ5xwQtLWuXNnj0855RSPY/4/Vl+HDh087tu3r8cvvfRS0q/QZb51T5tddtklaTvssMM8jud2nc8ojp///Oce6/w1S8dN95XjGlkzep9xxBFHJG16DtT5d+WVVyb9irmvl+5ZZGb2j3/8w+O4x6DuXfPyyy97XKrfBS2VXNtS6jo/9B71vPPOS/rpfjezZs3y+Mgjj0z6vfrqqx6X6udaivR3iY6pjpuZ2Q477ODx7373O4/1d2p8v0h/+82fP9/j999/P+mn5+LGjRsnbW+88YbHL7zwgseTJk1K+rGPzcp0fPUz1rE1M+vdu7fHcd/F1q1be6zjNn78+KTf1KlTPV62bFktj3jNsNIGAAAAAAAgg3hoAwAAAAAAkEF1nh6lS0hjWdIWLVp4HFNhdFlqr169PB4xYkTST1NrPvroI48/+eSTpN/mm2/ucdOmTZO2pUuXepxvSeO8efM8HjJkSNL29ttv5/x35SR+Ppq6Nnv27KRNx/u9997L+R5KlxGOHDkyadPSsjE97Z577sl32MiQvffe2+NWrVolbbq0nPKZhaOlYHUJ/ZZbbpn0y1U6Nc5ZTX2M59P999/f41133dXjfv36Jf3++9//1uTQK1pcpn3IIYd4rNfMe++9N+n3+eefF+04unfvnrS1bNnSY73eozh0ebiZ2QUXXOCxphebpfP2r3/9a7X/Hblpqd6YeqbpRkOHDvV44cKFRT0mHf/BgwcnbZ06dfI4pmscd9xxHuv2A6WqEGlnm266qcd/+MMfPG7evHnST89/mh561113Jf30e/D8888nbcy5NaNjENMC27Zt6/Gbb77pcUyb07axY8d6HFO89fdinM9/+9vfPH7kkUc8jukzX3zxhceLFi1K2r755ptq/x33vN/R+wjdRsUsnXPt2rXzOF4XVbyP0td6D7zzzjsn/f797397HH+P6tYsOtaFnuestAEAAAAAAMggHtoAAAAAAABkEA9tAAAAAAAAMqjB6uRbNWjQYI2TszTPTPMOzcz22msvj3V/G7O0XJeWBY55xZrbqHlquueCmVmjRo08jrmwmrPYrFkzj2OZNs1DvPjii5O2M8880wqpqqoqd8251VCIMSw2/Y5oGb6LLrooZ7/p06cnbbqPRqH3c1gDM6uqqroV4o1KYRxziWUXtQyw7odhZta1a1ePX3/99eIeWA2V4lzU852Z2YwZMzzecccd9ZiSfvn2sckl5mHrPNX31z3HzMy6dfthatRBOfCSnIt6PTIze+KJJzzWz7lHjx5JP82hL0SOtZZEfeWVV5I23Vsn5p/rXmaFUIpzsdBiKWctGatlkM3SkrGdO3f2uNDjspoyOxfjXhla8nXddddN2nQfDN37oBh7l+h5VMdRz+tm6X1v3759kza97hZCKc7FuO/Ftdde6/GwYcM8jntDqXzXSJ1vcb8h/S7pb4li0P/PeIzhdWbnYrw3Of744z0+44wzkrZ33nnH49/85jce59trVMc4jrf+Rsw3Vrn2AKzudTGV4lyMtt12W4///Oc/e6z3iWYrX+O+p/vKmJm99NJLHsffIFtssYXHeq8cx0z3ZY37gD311FMen3TSSR5rGfjq3jOPauciK20AAAAAAAAyiIc2AAAAAAAAGVTn9Th16XxMd9Dya3EJUa5Up7jsX5dmaxqVlrg1S5d3a6nGeIy6XPKAAw7IeUyxvDVqT1MAzj//fI9j+Vhd/nbUUUclbcuXLy/S0WFN6VLE+FqXE5utvLQQNafnJ10ibJamneUrjbhixQqPP/vsM4+ffvrppN+ECRM8jstSzznnnGr/bkw9ePDBBz3u1atX0hZTqSqJjs9ZZ52VtGl68O9//3uPi11muE2bNh7HlC29fsbSpiiMfGWe49JvpSkZH374YeEPrAzoeXPvvfdO2jp06OBxLKH98ssve1zsVAgd4+uvv95jvf81M7v99ts9LnQ6VDnQ3whmZnvssYfHmgYTS0Xrb5emTZt6vMkmmyT99HfGxIkTk7aDDz7YY01zLcZ3pxxKR8c0UN0uIW5/oSm7H3zwgcf5Upv0M9I0mNVBGffVo6lNp512WtI2YsQIjzUVNY7hG2+84fFNN93k8bhx45J+OofXWWedpE3TyXVLjh122CHpt9FGG3kc57re2+p5Rb9/hcBKGwAAAAAAgAzioQ0AAAAAAEAG1Xl6VD75lvDpsjNd5hSXLS5evNjjefPmeaw7R6+KLo/VJa/77rtv0k9TcO67774avz9Scad2XdKry+Li0sObb77Z40IvQUPx7LnnnslrXeo/c+bMpI00t9rTJaCnn3560hZTDb+n6VBmZrNmzfJ46NCh1f53s/zLjrWajaaibr755km/du3aeXzdddclbQMHDvS4HJZ6rw6tsnj44YcnbQsWLPBYz4fFoNfFMWPGeBy/S1pVIVZmRGHoNbN///5Jm45TnJfjx4/3OM51fEc/v/bt2ydtmpakS+XN0ip8mupZ05SJWB1HX8fl/LfddpvHO+20k8d6/2uWVtjBd/RzjXNHK8TOmTPH40GDBiX9Xn31VY9bt27t8f/+97+kn6YAx3QKTRvu16+fx8WuJFVKdKxGjx6dtOl1J6YDH3vssR7XtHpsvspPKIyYuqvp3iNHjkza9Bqn1cDiPZCm6tf0mpavYvSLL77ocZcuXXIeU9xWQF9rel2+iqy1wUobAAAAAACADOKhDQAAAAAAQAbx0AYAAAAAACCDMrWnTVZo3t2vfvWrnP00r/jTTz8t6jGVs2OOOSZ53bdv32r7PfPMM8nrM844w+OYy6j5rtpGrmr90FzQAQMGJG2aw33ppZfmbMPq2XjjjT2O+fRK94jRPbzM0vPf22+/Xe2/WRXdb+rOO+/0+OSTT0766Xl3t912S9q0dGrct6HcxBzoa6+91uNY0ve8887zuLZlSmuqefPmHvfs2dPjeE695557PK60/YfqSosWLTyO+66oWNY7XkOxMv3O6nfeLN23IJaLvuWWWzweO3asx3Gfk80228xjPS/Hffl0jGM53JYtW3qs8++SSy5J+n3xxReGlI7h4MGDkza93zj//PM9jtdF/cz1urhs2bKkn+5pE+n+Oahes2bNPI6/C/S+/vLLL0/adA+UmuK3QfHp/ltmZsOHD/c47nej+87o74J4DdM5q/dOca+9pk2bety5c+ek7ZBDDvG4T58+Hut5Nr5//L5MnTrVY91Pt9DfK1baAAAAAAAAZBAPbQAAAAAAADKI9ChbeTn6xRdf7LEuPV60aFHS78wzzyzugZWxDTfc0GNd4m+WptJo+diYVvPJJ594HJeg6RJnlj3WP11K3qNHj6RNlxQ/+eSTdXZM5W7LLbf0OJYn1Pnx/vvve3zggQcm/ebOnVvtv1kdOv+0LPURRxyR9NN0rvXWWy9p0+Ws//nPf2p1HKUiLhPWUsIx3eHuu+8u2nHE66KWvdWl/XqONjM7++yzi3ZMlUzHY++99/a4UaNGST+db7feemvStnz58iIdXXm67LLLkteHHXaYxzF1Su9pxowZ43G8/9Bx1BQPvZ8xS8/ZmioV30Ovn1dccUU1/xfIJV5ntLR6p06dPM5Xtrdhw4Y531+vmfH6qWOv53zm6A805URTpM3SdODp06ev8d/id0Jx6NyJ82jdddfN2aav9V4xjpPe52633XYexy0BdL7FNj1361yM9836t+PzgFGjRnkcS4oXEittAAAAAAAAMoiHNgAAAAAAABlUselRuvSqdevWSduwYcOq7Td69Oik35IlS4p0dOVJd/O+7rrrPI477OsytnPPPdfjd999N+mXbzkjSx3rn86dQw891OOY/jFz5kyPP//88+IfWIXQZaOxCpcuLdYqbLrrvVnulKh8y8Xz0aX88fyp6VFxWaqeo8s9PSqeD9daay2PP/roo6QtVisppFh9YeDAgR7r92Ly5MlJv7hsGIWhaVBHHnmkx3Eu6tLsm266qejHVc40PdTMrHv37h7rGJiZbbPNNh536NDBY025MUvPbXfddZfHWlnPzKx///4ex+pReg29+uqrPSa1ZtX0WvXee+8lbdtvv73Hp59+usex0tMLL7zgsVY6jNetpUuXVvt3zdL0q1133dXjBx54IO/xVxLdniLeN+p5b8iQIUnbnDlzPNZ7ynz3KXpNy5eqE9+D3xr56eej88bM7LnnnvN42223Tdp0PDTtKaaA6nYaep8bU5S0X77vkh5vvP/V4823XUcxsdIGAAAAAAAgg3hoAwAAAAAAkEE8tAEAAAAAAMigit3TRvPbNHfVLN1D4IMPPvB4woQJxT+wMqblhDXWsTAze/PNNz1+8MEHPSZ3tLTouGrJvphreuWVV3pc27LSWDmfvl+/fjnbFi9e7PEjjzzicdz7JpeY863vH8dQ90YZOnSoxy1btsz5/nGuV9L+YfH/XfepiOPTsWNHjzXfuqbjGOmc7dOnT9Kme3boeMdyq8zh4tDyxDH/X+lcefvtt4t5SBVHP89zzjknadNzos6PfHtgaBzPqXofetJJJyVtutfO2LFja3Lo+H96fpo4cWLSptdMLUccfyPo+fXLL7/0WM/BZmZffPGFxzvuuGPSpiXjr7nmGo+32GKLnMdbCXTuxH3VlM6Po48+Omk74IADPNbPL54PdT8o3QdnwYIFSb+4p5HSfeV0DyPdNxDfifdxe+yxh8fxmvaLX/zC465du3rcqlWrpN+nn37qsd7Xxr3E9G/Fc63SOXvmmWcmbddff73Hn332Wc73KCZW2gAAAAAAAGQQD20AAAAAAAAyqGLTo7Ss6n777Ze0aYk4LVVc2yXnlapp06bJ6/Hjx3usSxt1OZqZ2Y033uhxLLmp8i1xI5Wq/v3sZz/zeKuttvI4lgR+6KGH6uyYyllMM+zWrZvHscShvq7NMt58cy8uS+3bt6/HJ598cs5++p66zNXM7LHHHlvtYyxV8f9dS5z36NEjabvooos8fuuttzzWa1h8rUvE43dGr4uxpHGTJk081vMry8Drxj777OOxpkrFa919993ncby2onh0HGpzrxhTWC+88EKP9X7JzGzSpEkea0oGVs/jjz+evJ4yZYrH+rsgpunoefPdd9/1ePDgwUk/LQN8wQUXJG3Dhw/3eJNNNvE4pn9UWoqjpjM988wzHu+yyy5JP50vmspmZtamTZtq+2255ZZJPy21rtfIDz/8MOmn7x//ls4/TY87/PDDk36axoPvaIqR3udU9zoXnYv6m/O6665L+uW7Z9Xj2HfffXMeQxZSFVlpAwAAAAAAkEE8tAEAAAAAAMggHtoAAAAAAABkUMXsaRNz97XUW7NmzZK2qVOnevzkk08W98DKjOYNnnvuuUmb7nGi+d+PPvpo0m/cuHEe59svQf9Wfe5hk5XjqG8xZ/Sggw7yWPOAp02blvQjJ78wYlnKDTbYwOM4NrqnTdxLIZd833N9v969eydtukeV7osSaenUyy67LGmrpO+Ifg5m6d4HuoeNWbpvkcZxDwbdY0P3OXn11VeTfu+//77H+v2JNLc75v+jMOIYjhw50mO9n4n7F/3xj3/0OAs5+MhNz6m//OUvkzbdvyruTXTrrbd6zF6LtRevK7qPl5aR7tKlS9JP55ze5y5cuDDn3xo1alTyWksat23b1mMt/22WliGvtLEeMGCAxzfccEPSpufHdu3aJW0bbbSRx/n2/9LzqF7v9N/EfxfvpfTeVstK63nYzOzUU0/1mH3gCkevcVtvvbXHOhZR/PxPOOEEj3Wfqyz+nmOlDQAAAAAAQAbx0AYAAAAAACCDyjo9Spf99+rVK2nbf//9PY7L0bU0X6UtR1xTumRx0KBBSZsuK9Tyd5oOZbbyeHwvpnEUe+maHq/+7bhsXZfOPv/880lbJX1/YllpLXmoy1BfeOGFpB9L+AtjxYoVOdviXFm+fLnHmlYVv9s6B3QZsJZwNzM75phjPD7wwAOTNk0/1ff7+uuvk37333+/xzE9KovLVOvKvHnzPB46dGjSpqVNtbR6LBu7YMECj3UJfyxpqee5CRMmJG177bWXxzp2lVaStq7EOaYpADoftMysWVqCGNmm18Vhw4YlbXpdvOOOO5K2OXPmeFzJ58Y1FT+7xYsXe3z55Zd7HO9R9HVNP/+YknH++ed7fNVVV3m83XbbJf3at2/v8SuvvFKjv1UuNPVWU6XM0jGIvw10KwZN0z/xxBOTfloCXN8jjndMT1S5Us01XdnMrGHDhh6THlU4Ol8mTZrkcb4Ut7jtif67rJ9PWWkDAAAAAACQQTy0AQAAAAAAyKCyS4/S5fe6g/gtt9yS9NM0gKeeeippi+kbqDlNochXfUSXoMUUIk3X0Pdo3bp10k/TBuKu/bnSkuLO77qErmvXrknbSSed5PGmm27qcayAo6leuszVbOXvXTmLFdp0Wa9+7jNnzkz6kR5VGDHdSOdHTLXQsdIlw0888UTSb5999vG4f//+Hjdt2jTn+8U5pq91rF9++eWk3/HHH+9xrhTJSheXVWv1J43jGKh8VcB0efc///nPpG3XXXf1WL9rjFXh6Nho9SCzdBm+jttf/vKXpB9L70uH3i/FtJj58+d7HCvR5EuFzSWmkOj3KV/6RyXReVXoeRTPtX//+989HjhwoMeahmqW3lPGCmOVNG7x/kbF+/25c+d6fOWVV3o8e/bspJ+OgaYvxc91yZIlHscUcrXWWmt5vGzZsqQt3/Gj5vS3mJnZfffd5/GGG26Y89/pGJ5++ulJm24XkHWstAEAAAAAAMggHtoAAAAAAABkEA9tAAAAAAAAMqjs9rTRfMOxY8d6rPvbmKWl/Q455JCkjT02ai9f2Txt09zDm2++Oen34osverzzzjt7HEu4ff755x7PmjUradP9HfKVZ+zZs6fHugdLPN5877Fo0SKPY0neStrTJuaa6nhpzvHTTz9dZ8dUSWLO9Pjx4z3W8qVm6b5MJ5xwQrWxmVmjRo08jnsW5RLnh+6/oPs0HHfccUm/uC8Vai9f2cp8bXqei/tm6H4rmvNfSfsq1KVevXolr/Xz17l+//33J/2yXrIUPxg0aJDHzZo1S9p0/41PPvlkjf9W3NOGPTYKI98eYfnoeXP06NEe77LLLkm/Tp06edynT5+k7aGHHqrV364kek2L+5XqNa6me9XoPVH8d7q/27/+9a+kH3uN1Z5+/meddVbSpuXdVdxTaNiwYR7HfWxLCSttAAAAAAAAMoiHNgAAAAAAABlUdulR++67r8e69DQuDT3ttNM8/uCDD4p/YBVi6dKlHsclaLrcW8cjptXo6zhuSkv0xfKousRUlwHHJcG6tDWWrtXljFrW+6WXXkr6abpPLM1ZSc4+++zktS4b1e+Flt5D4cTl0ZqaF1ORtt9+e4+17GyUr3R0rr8d55GW9h4+fLjHMU2O5d3Z0qFDh5xtmppam/LDWDVNizDLnZ6m1yZkn5baPuWUUzxeZ511kn6NGzf2eIMNNkjadOl/vnR+vX+KZZE53xZGbT9H/Xd6T3nPPfck/Q444ACPDzrooKRtxowZHsfrLlam22KYmX388cce6+8O/W1hls7NON56b3vXXXd5fOONNyb92Haj9lq3bu3xYYcdlrTlOsfp9ihmZvfee6/HpTwWrLQBAAAAAADIIB7aAAAAAAAAZBAPbQAAAAAAADKo5Pe0ad68efL62muv9VhL1Gr5RLOVy0yjMHTPmP322y9p031njjnmGI+32mqrpN/666/vsY6vln0zS/dViHvVaBlu3Zsm7l904YUXehxLDr/zzjvVvkfMh6xtycdyoPn5u+++e9Km+aXXXHONx5QarRs6P+KeNlOnTvVY90uIe0jp91njuI/JvHnzPH744YeTNi1nqv0qba6UAt2Hqnfv3kmbnuc0j5/5XDj6Gce9pvR8+tprr3nM519adH8M3Zsoatu2rcc33HBD0nbyySd7PHfuXI9jWWG9V+F8m106brfeemvS1r59e49btGiRtOk9l17T4/5F+E6cb1oGWvdA2WKLLZJ+OmfjHmJXXXWVx1dffbXHeo3E6tN70TPOOMPjJk2aJP30vDZ9+nSPL7300qRfKe9jo1hpAwAAAAAAkEE8tAEAAAAAAMigkkyP0iXE48aNS9o23nhjj3WJYCxHXC5LpbJMy1KamU2bNq3auBhylSouxhLhSl52rEsYYwqiLtvWFJlK/rzqy8yZM5PXWvJ78ODBHmt5UbM0/W3WrFke6zJUM7MpU6Z4HNMMKQldOnTpcUzP0eXeugyc8S2O5cuXJ6/1fmb+/Pl1fTgoEE1b1WtkTBPXssNdu3ZN2gYMGOCxlhbWEsZmudNbkS06Ns8991zS9vjjj3vcs2fPpE1TWPWaHMt/83unevqZ9e/f3+M999wz6detWzePJ0+enLQ99thjHpOWVjh77LGHxwcffLDH8bedpryNGjWq2v9eTlhpAwAAAAAAkEE8tAEAAAAAAMigBquzZLJBgwaZWF/Zrl07j59//vmkrXHjxh4vWbLE45YtWyb9tK0UVFVVVZ/vs5qyMoYVamZVVVW3VXdbNcax/jAXywJzsRpacbFfv35JW5cuXTy+/vrrPdZKfXWtnOei3ueYmY0YMcJjrcb34osvJv1KMA2mYudix44dPZ44cWLSpulSmiJjZnbqqad6rFUu6zNVsZznYl2K6R/62+XOO+9M2t566y2P9TyglZDMVuucULFzMR+9LsZUsyyeb0txLmo6qJnZs88+67GeJ+Pnfdttt3k8dOhQj8sgJbDauchKGwAAAAAAgAzioQ0AAAAAAEAG8dAGAAAAAAAgg0pmTxstLfzoo4963KdPn5z/5o033vB4m222SdpiOc2sK8UcRayEfOEywFwsC8zFMsBcLAvMxTLAXCwO3etj2rRpSZuWpta9b15//fXa/jnmYhkoxbkY97SZM2eOx61atfL4vffeS/rpfjfLli0r0tHVC/a0AQAAAAAAKBU8tAEAAAAAAMign9T3AdSUpke1bt06Z78vv/zS46OPPrra/w4AAAAAdU1Le+fbpkJ/u/Tu3cMViOsAAAEvSURBVDtp099Fa5ASBdS7b775JnmtJewbN27s8ZAhQ5J+ZZYStUqstAEAAAAAAMggHtoAAAAAAABkEA9tAAAAAAAAMqhk9rRZsWKFx23atPF47bXXTvp99dVXdXVIAAAAAFBj+faxqalvv/22AEcC1D/9jW9mtv/++9fTkWQbK20AAAAAAAAyiIc2AAAAAAAAGbS66VGLzOydYhxIbVVIOlTuGuerL3NjWEEYx9LHGJYHxrH0MYblgXEsfYxheWAcSx9jWB6qHccGhcirBAAAAAAAQGGRHgUAAAAAAJBBPLQBAAAAAADIIB7aAAAAAAAAZBAPbQAAAAAAADKIhzYAAAAAAAAZxEMbAAAAAACADOKhDQAAAAAAQAbx0AYAAAAAACCDeGgDAAAAAACQQf8HfI5aWx9qkz0AAAAASUVORK5CYII=\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: 1.85 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": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0024, PSNR (dB) = 26.18\n",
      "MSE = 0.0032, PSNR (dB) = 24.93\n",
      "MSE = 0.0007, PSNR (dB) = 31.83\n",
      "MSE = 0.0033, PSNR (dB) = 24.81\n",
      "MSE = 0.0032, PSNR (dB) = 24.94\n",
      "MSE = 0.0006, PSNR (dB) = 32.45\n",
      "MSE = 0.0050, PSNR (dB) = 22.99\n",
      "MSE = 0.0053, PSNR (dB) = 22.74\n",
      "MSE = 0.0052, PSNR (dB) = 22.81\n",
      "MSE = 0.0040, PSNR (dB) = 23.94\n",
      "Wall time: 8 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.9772\n",
      "0.9769\n",
      "0.9898\n",
      "0.9713\n",
      "0.9576\n",
      "0.9923\n",
      "0.9354\n",
      "0.9535\n",
      "0.9424\n",
      "0.9638\n",
      "Wall time: 136 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": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9960\n",
      "0.9969\n",
      "0.9986\n",
      "0.9953\n",
      "0.9930\n",
      "0.9988\n",
      "0.9927\n",
      "0.9928\n",
      "0.9871\n",
      "0.9901\n",
      "Wall time: 140 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/tensorflow/models/blob/master/research/compression/image_encoder/msssim.py\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from scipy.ndimage.filters import convolve\n",
    "\n",
    "def _FSpecialGauss(size, sigma):\n",
    "  \"\"\"Function to mimic the 'fspecial' gaussian MATLAB function.\"\"\"\n",
    "  radius = size // 2\n",
    "  offset = 0.0\n",
    "  start, stop = -radius, radius + 1\n",
    "  if size % 2 == 0:\n",
    "    offset = 0.5\n",
    "    stop -= 1\n",
    "  x, y = np.mgrid[offset + start:stop, offset + start:stop]\n",
    "  assert len(x) == size\n",
    "  g = np.exp(-((x**2 + y**2)/(2.0 * sigma**2)))\n",
    "  return g / g.sum()\n",
    "\n",
    "\n",
    "def _SSIMForMultiScale(img1, img2, max_val=255, filter_size=11,\n",
    "                       filter_sigma=1.5, k1=0.01, k2=0.03):\n",
    "  \"\"\"Return the Structural Similarity Map between `img1` and `img2`.\n",
    "  This function attempts to match the functionality of ssim_index_new.m by\n",
    "  Zhou Wang: http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "  Returns:\n",
    "    Pair containing the mean SSIM and contrast sensitivity between `img1` and\n",
    "    `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  img1 = img1.astype(np.float64)\n",
    "  img2 = img2.astype(np.float64)\n",
    "  _, height, width, _ = img1.shape\n",
    "\n",
    "  # Filter size can't be larger than height or width of images.\n",
    "  size = min(filter_size, height, width)\n",
    "\n",
    "  # Scale down sigma if a smaller filter size is used.\n",
    "  sigma = size * filter_sigma / filter_size if filter_size else 0\n",
    "\n",
    "  if filter_size:\n",
    "    window = np.reshape(_FSpecialGauss(size, sigma), (1, size, size, 1))\n",
    "    mu1 = signal.fftconvolve(img1, window, mode='valid')\n",
    "    mu2 = signal.fftconvolve(img2, window, mode='valid')\n",
    "    sigma11 = signal.fftconvolve(img1 * img1, window, mode='valid')\n",
    "    sigma22 = signal.fftconvolve(img2 * img2, window, mode='valid')\n",
    "    sigma12 = signal.fftconvolve(img1 * img2, window, mode='valid')\n",
    "  else:\n",
    "    # Empty blur kernel so no need to convolve.\n",
    "    mu1, mu2 = img1, img2\n",
    "    sigma11 = img1 * img1\n",
    "    sigma22 = img2 * img2\n",
    "    sigma12 = img1 * img2\n",
    "\n",
    "  mu11 = mu1 * mu1\n",
    "  mu22 = mu2 * mu2\n",
    "  mu12 = mu1 * mu2\n",
    "  sigma11 -= mu11\n",
    "  sigma22 -= mu22\n",
    "  sigma12 -= mu12\n",
    "\n",
    "  # Calculate intermediate values used by both ssim and cs_map.\n",
    "  c1 = (k1 * max_val) ** 2\n",
    "  c2 = (k2 * max_val) ** 2\n",
    "  v1 = 2.0 * sigma12 + c2\n",
    "  v2 = sigma11 + sigma22 + c2\n",
    "  ssim = np.mean((((2.0 * mu12 + c1) * v1) / ((mu11 + mu22 + c1) * v2)))\n",
    "  cs = np.mean(v1 / v2)\n",
    "  return ssim, cs\n",
    "\n",
    "\n",
    "def MultiScaleSSIM(img1, img2, max_val=255, filter_size=11, filter_sigma=1.5,\n",
    "                   k1=0.01, k2=0.03, weights=None):\n",
    "  \"\"\"Return the MS-SSIM score between `img1` and `img2`.\n",
    "  This function implements Multi-Scale Structural Similarity (MS-SSIM) Image\n",
    "  Quality Assessment according to Zhou Wang's paper, \"Multi-scale structural\n",
    "  similarity for image quality assessment\" (2003).\n",
    "  Link: https://ece.uwaterloo.ca/~z70wang/publications/msssim.pdf\n",
    "  Author's MATLAB implementation:\n",
    "  http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "    weights: List of weights for each level; if none, use five levels and the\n",
    "      weights from the original paper.\n",
    "  Returns:\n",
    "    MS-SSIM score between `img1` and `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  # Note: default weights don't sum to 1.0 but do match the paper / matlab code.\n",
    "  weights = np.array(weights if weights else\n",
    "                     [0.0448, 0.2856, 0.3001, 0.2363, 0.1333])\n",
    "  levels = weights.size\n",
    "  downsample_filter = np.ones((1, 2, 2, 1)) / 4.0\n",
    "  im1, im2 = [x.astype(np.float64) for x in [img1, img2]]\n",
    "  mssim = np.array([])\n",
    "  mcs = np.array([])\n",
    "  for _ in range(levels):\n",
    "    ssim, cs = _SSIMForMultiScale(\n",
    "        im1, im2, max_val=max_val, filter_size=filter_size,\n",
    "        filter_sigma=filter_sigma, k1=k1, k2=k2)\n",
    "    mssim = np.append(mssim, ssim)\n",
    "    mcs = np.append(mcs, cs)\n",
    "    filtered = [convolve(im, downsample_filter, mode='reflect')\n",
    "                for im in [im1, im2]]\n",
    "    im1, im2 = [x[:, ::2, ::2, :] for x in filtered]\n",
    "  return (np.prod(mcs[0:levels-1] ** weights[0:levels-1]) *\n",
    "          (mssim[levels-1] ** weights[levels-1]))\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = MultiScaleSSIM(original, contrast, max_val=1.0)\n",
    "    print(f'{a:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Deep Autoencoders in Keras - 784-256-128-64-128-256-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
}