Simple Autoencoders in Keras - 784-64-784-tf2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install 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": [
      "1595257645\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 50s - loss: 0.6931WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.195999). Check your callbacks.\n",
      "235/235 [==============================] - 11s 48ms/step - loss: 0.2452 - val_loss: 0.1621\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 13s 57ms/step - loss: 0.1432 - val_loss: 0.1262\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 14s 60ms/step - loss: 0.1178 - val_loss: 0.1081\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 20s 83ms/step - loss: 0.1035 - val_loss: 0.0969\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 19s 81ms/step - loss: 0.0943 - val_loss: 0.0897\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 15s 64ms/step - loss: 0.0880 - val_loss: 0.0846\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 17s 72ms/step - loss: 0.0837 - val_loss: 0.0811\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 17s 71ms/step - loss: 0.0808 - val_loss: 0.0788\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 15s 66ms/step - loss: 0.0788 - val_loss: 0.0773\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 18s 77ms/step - loss: 0.0775 - val_loss: 0.0762\n",
      "Epoch 11/50\n",
      "235/235 [==============================] - 19s 80ms/step - loss: 0.0765 - val_loss: 0.0753\n",
      "Epoch 12/50\n",
      "235/235 [==============================] - 16s 70ms/step - loss: 0.0758 - val_loss: 0.0747\n",
      "Epoch 13/50\n",
      "235/235 [==============================] - 35s 149ms/step - loss: 0.0743 - val_loss: 0.0737\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 31s 132ms/step - loss: 0.0741 - val_loss: 0.0734\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 23s 99ms/step - loss: 0.0740 - val_loss: 0.0733\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 21s 89ms/step - loss: 0.0739 - val_loss: 0.0732\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 25s 106ms/step - loss: 0.0737 - val_loss: 0.0731\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 26s 112ms/step - loss: 0.0736 - val_loss: 0.0731\n",
      "Epoch 22/50\n",
      "235/235 [==============================] - 27s 114ms/step - loss: 0.0735 - val_loss: 0.0730\n",
      "Epoch 23/50\n",
      "235/235 [==============================] - 30s 126ms/step - loss: 0.0735 - val_loss: 0.0728\n",
      "Epoch 24/50\n",
      "235/235 [==============================] - 28s 117ms/step - loss: 0.0734 - val_loss: 0.0728\n",
      "Epoch 25/50\n",
      "235/235 [==============================] - 29s 123ms/step - loss: 0.0733 - val_loss: 0.0728\n",
      "Epoch 26/50\n",
      "235/235 [==============================] - 31s 130ms/step - loss: 0.0733 - val_loss: 0.0727\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 32s 136ms/step - loss: 0.0732 - val_loss: 0.0727\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 31s 133ms/step - loss: 0.0732 - val_loss: 0.0726\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 32s 138ms/step - loss: 0.0732 - val_loss: 0.0726\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 32s 136ms/step - loss: 0.0731 - val_loss: 0.0725\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 32s 136ms/step - loss: 0.0731 - val_loss: 0.0726\n",
      "Epoch 32/50\n",
      "235/235 [==============================] - 37s 156ms/step - loss: 0.0730 - val_loss: 0.0724\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 37s 158ms/step - loss: 0.0729 - val_loss: 0.0725\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 37s 156ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 38s 162ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 36s 154ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 35s 151ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 34s 146ms/step - loss: 0.0729 - val_loss: 0.0724\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 37s 157ms/step - loss: 0.0728 - val_loss: 0.0724\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 36s 151ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 34s 144ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 38s 164ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 38s 161ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 37s 157ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 37s 156ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 40s 169ms/step - loss: 0.0727 - val_loss: 0.0723\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 38s 160ms/step - loss: 0.0727 - val_loss: 0.0723\n",
      "1595257645\n",
      "Wall time: 23min 24s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import os\n",
    "import time\n",
    "\n",
    "CURRENT_TIME = int(time.time())\n",
    "print(CURRENT_TIME)\n",
    "\n",
    "try:\n",
    "    os.makedirs(f'autoencoder_mnist/logs-{CURRENT_TIME}')\n",
    "except Exception:\n",
    "    pass\n",
    "\n",
    "# Run this in a terminal\n",
    "#tensorboard --logdir=autoencoder_mnist/logs\n",
    "\n",
    "CALLBACKS=[\n",
    "    TensorBoard(log_dir=f'autoencoder_mnist/logs-{CURRENT_TIME}')\n",
    "]\n",
    "\n",
    "# Fit data\n",
    "autoencoder.fit(x_train, x_train, epochs=50, batch_size=256, shuffle=True, validation_data=(x_test, x_test), callbacks=CALLBACKS)\n",
    "\n",
    "print(CURRENT_TIME)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 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+WH4yJAAAURElEQVR4nO3bfXTP9f/H8dd85rPrmW1msjFMVxSFM+KoITpYKlO6ni505KgTlYij4+ToVGQdk1Nz6ii16rASxVoXJ2RCGFtliG3GbDJ2kdnF5/fH989fPr/X45xPv++rzv32993Tp70/76unFeTz+QwAAAAAAADc0uG//QEAAAAAAADwv7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwULASe71eX1hYmHXfrVs36cP88ccfUt+pUyfrtqGhQZrd3t4u9RcuXJD6ixcvSr3P5wuS/sBlREZG+mJjY637uro6aX5ISIj6kawlJSVJ/ZkzZ6RePebKz9EYY0pLS2t9Pl8X6Q9dRlhYmC8qKsq6b2xslObHxMRIvXLc29rapNmtra1Srx539efY3NwckHPR6/X6QkNDrXuPxyPN79BB28kHB9vfDnw+nzS7ublZ6oOCtB9x586dpf748eMBOxc7derkS0hIsO7V+1x8fLzU19fXW7fqfUu55xqjn+uXLl2ybgN5LoaGhkrX05aWFmm+egyVe5HyTGaM/ix07tw5qVefAWprawN2LkZERPiUe5d6HWtqapJ65Xw5efKkNFu55hijn4uRkZHWbU1Njblw4UJAzsWOHTtK90XlvDXGmIiICKmvrKy0bnv16iXN/vXXX6W+a9euUu/1eqW+vLw8YOdiZGSkLy4uzro/deqUND86OlrqlXNRfe9Rr8Hqu4NyHI8fP25qa2sDci6Gh4dL11P1GVV951KuC+p3X31Hqq2tlXp1H3Ls2LG/PBelpU1YWJgZPny4dT9//nxlvMnLy5P6cePGWbc7d+6UZqs358LCQqkvKSmR+kCJjY01c+bMse6/+OILaX5KSor4iey99tprUr9y5UqpV18sp0yZIvUDBw48If0BP6KiokxmZqZ1v3fvXmn+pEmTpL5Hjx7WrfrSoF4cc3JypP7mm2+2bgsKCqTZ/oSGhpq0tDTrXn04VW5wxmgvAcrLtTH/eZBQqAsn5Vwwxphp06YF7FxMSEgw2dnZ1v0HH3wgzX/00Uel/vvvv7duv/vuO2n22LFjpV59ECovL7dut27dKs32JyoqSrrmVVdXS/Mff/xxqVeWaQMHDpRmb9++Xeo/+eQTqe/Tp4/U5+bmBuxcjImJMTNmzLDu1eWbeh+dMGGCdbtgwQJp9vTp06VeXdAOGzbMun3hhRek2f6EhoZK3+lRo0ZJ84cOHSr1zz//vHX77rvvSrNHjBgh9VlZWVKfnJws9TNnzgzYuRgXF2defPFF6/7ll1+W5o8ePVrqlXNxw4YN0uwBAwZI/X333Sf1ynEcPHiwNNufmJgY89hjj1n36j/qqM8H/fr1s27Vf+TfvXu31Ofm5kr9vHnzpH7q1Kl/eS7yv0cBAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOChYiX0+n2ltbbXuc3NzpQ/j9XqlfseOHdbt559/Ls1eunSp1K9YsULq77jjDuv2+++/l2b7U1NTY1avXm3dT5kyRZpfVFQk9QMHDrRus7KypNnqMR83bpzU9+3bV+oDqb293TQ1NVn3yrlijDGdO3eW+mnTplm3Q4cOlWYfPnxY6ktLS6V+7ty51q3H45Fm+9Pc3GyOHDli3V9xxRXS/DNnzkj9jBkzrNvdu3dLs2+77Tapz8nJkfr8/HypD6Tm5mZz7Ngx676srEyav2XLFqlvaWmxbleuXCnNrqiokPorr7xS6g8ePGjd/vTTT9Jsf1pbW01dXZ11/+GHH0rzo6Ojpf7555+3btVrUq9evaT+q6++knr1XA+klpYWU1VVZd3HxcVJ84cMGSL1n376qXWrfqeSk5OlvrGxUeq3bt1q3V68eFGa7U9sbKx58MEHrXv1mfPbb7+V+sLCQus2Oztbmq0+Q2ZmZkr966+/LvWBFBYWZq6++mrrfsKECdL8kydPSn1xcbF1u2jRImm2+q6hPPcZY0x4eLh1q7yj/1/Onj1r1q1bZ92vWbNGmv/ss89K/YgRI6zb6upqafbIkSOlvlu3blKv3Av84TdtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcFKzEQUFBJjjY/o+kpaVJHyYvL0/qp06dat0OHjxYmp2cnCz1I0aMkPr4+HipD5SkpCSzfPly6760tFSaP378eKlXjnn//v2l2RMnTpT68vJyqd+0aZPUB1Jra6s5e/asdf/LL79I89XvZ2trq3VbVFQkzW5oaJD6rKwsqff5fNZtS0uLNNuf8PBwc8MNN1j3p06dkuYnJiZKfUFBgXUbHh4uzf7uu++kvnv37lLv8XikPpDi4uLM/fffb90fPnxYmh8VFSX1+fn51m1dXZ00+7333pP6jIwMqVe+sxcuXJBm+9Pa2mpqa2ute/WYpKSkSL1ybVeOtzHGjBo1SurfeustqQ8LC5P6QIqIiJCeO0tKSqT5J0+elPqEhATrdv/+/dLscePGSX12drbUf/vtt9ZtfX29NNufP//80+zbt8+6T01NleZ37NhR6qdNm2bdxsXFSbMPHjwo9R988IHUK89lgdbS0mKqq6ut+8zMTGn+2rVrpV55BlGfnQ4dOiT1c+fOlXrlO1hZWSnN9ickJMT06dPHul+8eLE0Xz13lXfX66+/XpqtvsNv2LBB6kNDQ6X+cvhNGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAABwUrsdfrNd27d7fuExISpA9z7bXXSv2WLVus2549e0qz169fL/X79u2T+rS0NOu2vb1dmu1PVVWVeemll6z7Ll26SPPT09Olvrq62rq94oorpNnPPPOM1H/++edSHxwsnT4BFR0dbcaNG2fd33777dL8oUOHSn1KSop1e+7cOWl2eXm51A8bNkzqi4uLrVuPxyPN9qepqUm6big/Y2OMue6666S+tLTUug0NDZVmDxo0SOoLCgqk/ujRo1IfSOo1tUePHtJ89fqfmppq3fbv31+afdddd0n9ww8/LPXKfbetrU2a7U9SUpJZunSpdX/NNddI86dOnSr1iYmJ1m1RUZE0+7fffpP6wYMHS31VVZXUB1J1dbVZvny5dT979mxpfn5+vtQrz8D19fXS7L59+0p9ZGSk1Cvnuvqd8sfr9Ur3uvfff1+ar94XleePe+65R5pdVlYm9Zs3b5b63r17S30gtbS0SM/36ndo9OjRUl9XV2fdqtf3iooKqb///vulXrk2BPJ9sb6+3nz99dfW/auvvirNP3DggNQrz9/qe8b06dOl3uv1Sv2QIUOk/nL4TRsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcFCwEsfHx5vHH3/cuv/oo4+kD9Pc3Cz1U6ZMsW4/++wzaXZoaKjUL168WOq/+eYb6zYoKEia7U9YWJjp16+fdT98+HBp/oIFC6R+9uzZ1m1aWpo0+9VXX5X68+fPS/3YsWOlPpAqKirMrFmzrPuYmBhpfkZGhtQXFhZatx9//LE0e9KkSVJfUlIi9SkpKdat1+uVZvvTqVMnM3HiROte/e86cOCA1D/00EPW7bp166TZe/bskfpt27ZJ/Zw5c6S+uLhY6v1JTEw08+bNs+7z8/Ol+RUVFVK/efNm67aqqkqafezYMak/evSo1Hfu3Nm6bWlpkWb7U1ZWZiZMmGDdK/dQY4w02xhjtmzZYt0WFRVJs++9916pV78jqampUh9I4eHhZtCgQdZ97969pfk7d+6U+uzsbOt29erV0uwzZ85I/aOPPir1p0+ftm5bW1ul2f40NTWZn3/++W/7u5cuXSr1q1atsm6V5yBjjMnLy5P6zMxMqX/nnXekPpDCw8PNjTfeaN2HhYVJ8w8dOiT1yvm1aNEiaXZERITUq5T7osfjCdjfO2jQILNr1y7rPjc3V5q/du1aqc/JybFuN27cKM1OT0+X+jFjxkj91q1bpf5y+E0bAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHBStxTU2NWb16tXUfEhIifZisrCyp/+2336zb9vZ2aXZ4eLjUFxcXS31kZKR16/F4pNn+BAcHm/j4eOt+wYIF0vzGxkap37t3r3UbExMjzU5LS5P633//Xeo3bdok9YF05ZVXSufiAw88IM1PTEyU+tTUVOv24YcflmZXVlZKfVJSktSXlZVZt83NzdJsf4KCgkzHjh2t+5aWloD93X/lxIkT1u25c+ek2cp/pzHG3HLLLVK/Zs0aqQ+k4uJi061bN+t+xowZ0vz169dL/XvvvWfd5uXlSbMPHz4s9QMHDpT6AQMGWLdHjx6VZvvTt29fs3btWuv+ueeek+b37NlT6rt3727dxsbGSrO3bdsm9cp32xhj2trapD6Q2tvbTUNDg3Wfn58vzV+4cKHUd+nSxbodM2aMNFs9t06fPi31yn3U6/VKs/1pb283Fy9etO7V94bhw4dLvXJvUZ4ljDGmvr5e6g8dOiT1R44ckfpAOnPmjFmxYoV1HxwsvY5K10hjjCkpKbFuT506Jc2OioqS+ujoaKlPSEiwbtVnLX/Ky8vNrFmzrHv1GKrPtOPHj7dut2/fLs1Wn7Pi4uKkPlDv8fymDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg4KV2Ov1mqSkJOu+qalJ+jBHjhyR+sOHD1u3RUVF0uxbbrlF6p944gmp//LLL61br9crzfanra3NnD9/3rpPT0+X5n/00UdS36GD/d6woKBAmt21a1epr62tlfqFCxdKfUZGhtT709zcbI4fP27dK983Y4y56aabpP7JJ5+0bletWiXN9vl8Up+amir1CQkJ1u2lS5ek2f506NDBhIaGWvfFxcXSfOU8N8aYhx56yLrdv3+/NHvmzJlS/80330h97969pf7UqVNS70+PHj3MCy+8YN0/9dRT0vwlS5ZIvfId/eGHH6TZV111ldSr39kdO3ZYtzU1NdJsfxobG83OnTut+2XLlknz7777bqm/9dZbrds9e/ZIs5WfsTHGvPbaa1JfUVEh9YEUGxtr7rvvPuv+7bfflua/8cYbUp+WlmbdnjhxQpqdmJgo9fn5+VKfl5dn3ar3aH/Onz9vNm7caN0PGDBAmj99+nSpz87Otm7VZ4/CwkKpV64LxujHfOTIkVLvj8fjMbGxsdZ9eXm5NH/y5MlS37FjR+u2qqpKmv3nn39KfWVlpdSvW7fOulXfu/1pa2sz9fX11r1yvTPGmJCQEKlXngs9Ho80W322GTt2rNTv3btX6i+H37QBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHBQsBKfO3fOfPrpp9b9m2++KX2YRx55ROrr6uqs27lz50qzg4OlH405ffq01BcWFlq39fX10mx/kpOTzbJly6z7a6+9VpqflJQk9SkpKdat8t0zxpiJEydK/axZs6T+xx9/lPpACg4ONjExMdb9K6+8Is0fNWqU1Pfp08e6nT9/vjS7qKhI6idNmiT12dnZ1m17e7s025/Gxkaze/du6/7uu++W5t9xxx1SHx4ebt12795dmq32GRkZUl9bWyv1O3bskHp/Ll26ZE6cOPG3/d2TJ0+W+hEjRli3+fn50ux3331X6vft2yf1d955p3X7/vvvS7P9CQkJMb1797bu16xZI82fPXu21AcFBVm3S5YskWZ/9tlnUn/jjTdK/SeffCL1gdTU1GT27t1r3cfHx0vzc3JypH7Xrl3WrfIsZMx/rjuK7du3S/2FCxes27a2Nmm2P/Hx8eauu+6y7isrK6X5ra2tUp+enm7dqt8P9Zg3NDRIfVlZmdQHUkhIiPRcmJubK81Xjosx2jOq+t6zfPlyqZ82bZrUK8+dgXy28Xg80ntG165dpflPP/201Kvvl4rExESpX7p0qdSHhIRI/eXwmzYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA4K8vl89nFQUI0x5sTf93FwGT19Pl+XQAziGP5XcRz/+TiG/w4cx38+juG/A8fxn49j+O/Acfzn4xj+O/zlcZSWNgAAAAAAAPj/wf8eBQAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOCg/wEXNWM/y2O7UAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de7xVYx7H8ScjFV2kojRJESEKTe6XIYMIhZcpYVxjNMLkbtxmjEHzch1KqJfGrXqNS7m9xrVyGakwmSiVJCmdlFEMY8784eXn+/yctc46u31O6+z9ef/12z3P2We3nv2stc56Pb/n16CysjIAAAAAAAAgX9Zb1x8AAAAAAAAAP8ZDGwAAAAAAgBzioQ0AAAAAAEAO8dAGAAAAAAAgh3hoAwAAAAAAkEM8tAEAAAAAAMih9WvSuUGDBtQHX0cqKysbFON9GMN1anllZWWbYrwR47juMBdLAnOxBDAXSwJzsQQwF0sCc7EEMBdLQpVzkZU2QN1ZuK4/AIAQAnMRyAvmIpAPzEUgH6qcizy0AQAAAAAAyCEe2gAAAAAAAOQQD20AAAAAAAByiIc2AAAAAAAAOVSj6lEAAAAAgGwaNIgL+lRWUpgHQM2w0gYAAAAAACCHeGgDAAAAAACQQ6RHoc40a9bM4iFDhlg8YMCAqF/79u0t/t///mdxo0aNon6zZs2yeNCgQVHbBx98UOV7oH7xS4q/x9Lidathw4bR6/XX/+FS8uWXX9b1xwEAIFf0/oV7lnVjvfV+WJtQjL8F9P2K9Z5AVqy0AQAAAAAAyCEe2gAAAAAAAOQQD20AAAAAAAByiD1tUFSa79mpU6eo7dlnn7W4Y8eOFiftW1KdXr16WTxp0qSorXfv3hZ//PHHBb1/Ep/Tqshvrbmf/OQnFm+33XZR29ChQy2eOHGixU888UTU79tvv62lT1fedGwuvvhii88999yo34IFCyzu06dP1FZRUWFxoXn97A1QPT0v+XOqnpc4fuVD528I8b5yq1evtvibb76ps88ElAvOteue7rf33//+N2orxv361ltvbbHu5/f5559H/b766iuL27RpE7Xp30qffvqpxQsXLoz6/ec//1m7D4t6j5U2AAAAAAAAOcRDGwAAAAAAgBwiPQpF1aNHD4t9Cstmm21msS7f12WDIcSpTi+++KLFxxxzTNRvjz32sLhz585R28CBAy0ePnx4lo+emV9yztLy7/iUjKxLg/V49u3bN2rbfffdLZ46darFPkWN9Kji8Mf16KOPtvh3v/udxY0aNYr6adrF3nvvHbU99thja/05VDkvOfdzbMMNN7S4S5cuFm+00UZRv7lz51qs6Wo1mTf6uzX2Y6Vt/v1JH619ej494ogjojadw/fcc4/Fd911V9SPa1r1GjduHL3u2bOnxa1bt7Z41qxZUb+PPvrI4rR0h6T5FkL2dEedm23bto3a9NyhKeSa4lHd+6Pw9H7FMS4uHZO0c5n20/PmxhtvHPVr1aqVxYceemjUdtFFF1ncsmVLi/11Uf/O8ddB/zfF92bPnh291vthrqXliZU2AAAAAAAAOcRDGwAAAAAAgBzKVXqULifzS8t0KRjLwvLDLw3VcdPlt96KFSss3nXXXaM23TFdl436ClEzZsywWNMzQoiXMxZaeSZpeTJLx6tW6BJf/c74cdTloS+88ILFvgoAimPLLbeMXt98880Wb7DBBpneo1u3btFrTZPMOnc4x2ejaVCaWuhTN0aNGmXx8uXLC/pdOr+1Ikf79u2jfrq0fM6cOVHbmjVrCvrd5SIt1SLr+bVp06YW6/wNIYSf/vSnFp9//vkW33fffVE/rnFV0zSG0047LWobNmyYxVo55uqrr476LVq0yOJCr5lJP+e/PzreF1xwQdSmc1PHnzn6Hb0v2XzzzaM2rU76/vvvWzxt2rSo39dff22xH7Ok1Bwv7Vqo71mMc0epyPr/1XuaCy+80GJNeQohTgf3f5umVW1UaWOs9D122GGHxM/rt5VAeWClDQAAAAAAQA7x0AYAAAAAACCHeGgDAAAAAACQQ3W+p03Dhg0tPvjgg6O2c845x+IWLVpEbW+//bbF48aNs3jBggVRvy+++KLKWHNLQ4jzRH2uob5OyxnVNl+6sVz2ZPC5o4sXL7ZYxyyEuAzmgAEDLP7ggw8y/a7u3btHr5s0aWKx3+NE99EoRt54ueUEFyJrye+0vPvddtstanvjjTcsXrJkSbXvjWw0D3uLLbawePz48VE/nbPKH389Z5577rlR22abbWbx8OHDLdbStyFQtr0QnTp1srhPnz4Wv/baa1G/ZcuWWZx2bcq6L4K+x0477RT1O/bYYy32+3nMmzevyvdD4fweC7rfht9vSPuuWrXKYvYxyaZr164WX3LJJVFb8+bNLX711Vctfv7556N+/l60mPxeVjfddJPFO+64Y9Sm5wRf5rtc6fw48sgjLR49enTUT/dr/PDDDy0+6qijon7vvPNO4u/S859e+/x8TvqZmrSVGz2GacdF9/DTfWx0XzDPXz917DTWfa1CCOG9996z2M+3TTfd1GLdp+7++++P+tWXfWyS7iOK8R1N20dVz8H77bdf1O+AAw6w2P8tqcf8rbfesvivf/1r1G/69OkWr169OmrT70VtzkVW2gAAAAAAAOQQD20AAAAAAAByqE7So3Q5ky5D8ulRusxalzmFEC9L1Z/zaTG6rE2XL61cuTKxn08B+OyzzyzWJW1t27aN+uny0t/+9rdRmy6LLKdli5rC4peKatm8pUuXZno//Zm//OUviW1+KeKsWbMyvT/WTlJZdE/ngE9H1LSObbbZJmobMWKExZT5Lpw/5rpUfuzYsRbrebaqn/teWvnSli1bRm2DBw+2WMtSa4nNEEKYMGGCxeWSXlpTfun8KaecYrGmuY0ZMybq59N3v+fnrI53Wrqajo9Pu9h9990t9qVytTwuakbHSuefppyHEMLpp59usZ+/Oqaajk5qYtX88RsyZIjFm2yySdS2cOFCi7X8t7/3VGmpMCrt2qptmpoYQgi/+MUvLF6xYkXU9vDDD1us19asac6lSP8WGDlypMV+uwbVoUMHiy+//PKo7bzzzrPY3/MmXeO49tWcpmCHEN+D6BYa/h5SX2uam6ZNhZC+7YNuxfDSSy9ZXFFREfXTFFR/vtU5ljW1K8+K/bm11Lm/R73xxhst3mOPPSzWvw/9Z9L38/R+RlMkQ4i38njqqaeitlGjRlmsqf/Fns+stAEAAAAAAMghHtoAAAAAAADkEA9tAAAAAAAAcqhO9rTRXDItdfbggw9G/TSXT0ughRDnA2rO7c477xz1a9asWZWfoV27dtFrzVH1eY5aKrxLly4W654BIYSw3XbbWax55CGEcP7551us+YvllC/sc6j1/5617OygQYMs9vsjaL+pU6dGbTqGqBv+u5z03fb5pEcffXTie86cObPa9wshfW8dcsR/fF7U/aH0PJa0h00I6cdfz6HffPNN1KZ7bmh599tvvz3qp2WqNb8cP9DrUQjxflBaDnTSpElRv6Q9S9LKZ6btc7L++j/cOvjSmloON20/D9SMzj8dN51TIYSw1157Jb6H7sOn5UvTlNM9i9epU6fotV6r/H2jlgDXPTA8PZ463/x1So+zb9P3aNWqlcW///3vo356Pr/vvvuitrlz51b5u8qJv97p/mt+z6Is79G/f/+oTc+NQ4cOjdoeeeQRi/01E9XbeuutLZ48eXLU9sYbb1h83HHHWeznke5Dtdtuu1ns71H12urHSq+TOo8KnVPlur+Y3lOEEEKvXr0svuaaayzWfWtCCKFJkyYW6zH3Jbm1lPfGG28ctbVv395iPSf7Pcd22GEHi7t16xa16V44xx9/fOLnWFustAEAAAAAAMghHtoAAAAAAADkUJ2kRyktPfqPf/wjaps2bZrFaakWunRNy+15q1atstiXhNOfmzFjRtSmqTWnnXaaxddff33UT5eoTpkyJWpLSskop2WoWf+vfvm1pq5deumlif001e6ss86K2igPXb2sJbqLzae5bb/99hZryfgQQvjkk08yvWdSOdxypsfELynt2bOnxWkpUUrn1Lx586I2nae65DiEOG1ASyj65edaunHgwIFRW7kuGQ4hhMaNG1t82WWXJbbdfPPNFvt5lMRfp3Tpd9o80lRhn6L8+eefW6wlMlEzacdfl22fdNJJUdtGG22U+HNanjapDHwI6deGUqfH1p+HdCm+fs9D+PE9YBId10LvUzTl9IorrrDYn1PffPNNi//whz8kfo5y5Ut5n3HGGVX20+0ZQojvS7S8tH8/3ebh2muvjdr+9a9/WTx79myL0651fl4m3feU4tj6FO9nnnnGYh2DEOJ5pdc4f1z0WLOlQt3TVOoHHnggajv44IMtTiuDrufhkSNHWnz11VdH/fTvRb1vCiGEXXbZxWK9N9a//0OIy4371Kl9993X4ubNm1vszx1rOzdZaQMAAAAAAJBDPLQBAAAAAADIoTpPj1J+GWDWJfC6zEl3wA8heenRp59+Gr3+5z//mfgzuuRQl175fhUVFRY///zzmT4HvqNLy5o2bRq1XX755RZrZQz//SgkHQA/KHZKUdbl/H7JuS579UvMddf+rL+bufcdTXu66qqrojZfHeF7Pl1m6dKlFmtFvMceeyzqp+Pkl3BfeeWVFmuFBt2xP4QQDj/8cItPPPHEqG3MmDEWl/r4+uOnqW0HHHBA1DZr1iyL77jjDouzVkyrSWUgre6g3yddChxCCO+8847Feq1G9bJ+t3X+9uvXL2rTMfXpN7p8POt3pNTnm6epR927d4/adA74NLStttrKYq2cmXb8kiqC+dc+NeSPf/yjxSeccILFfr7peTQtHa6c6HHVdN0Q4vQyvafU61YI8d8Thx12mMUPPfRQ1E+/Sz4tXH/uvffes7jQVOBSnKc6ViNGjIja9LozZ86cqE23S/j666+rfD8va3XSrD9XzlX3kjRq1Ch6PX78eIsPPfTQxJ/Tv/PPPvvsqO2VV16xWM9/acfbVwCbP3++xfq98lWsfYUrpc8K0p4brC1W2gAAAAAAAOQQD20AAAAAAAByiIc2AAAAAAAAObRO97Qphqz5Ylnzt0OIy4ENGzYs8T3uvPNOi1euXJn5/RHn5PscxcGDB1usx/y5556L+l1zzTUW+zzgpBzUcs8rLUaJyELKwWou69FHHx21aX7p2LFjo7as+d3lPq5V0fKjWlbd0znmyy7q3Pz3v/9tcdZ9GkKI85Hvuecei335ai2nq+W/Q4hLFS9btizxd5cCX45S9wTybbfccovFafs/6f5Gmpet+f6en+fbbrutxbr/kB/v0aNHW1xoSeNyUsjeYq1atbK4Y8eOif0+/PDD6LWWGU5TzudTPR/6PROV34tvwoQJFo8bN87iqVOnRv10vxsd+w4dOkT9dD+/M888M7FNx+r++++P+i1YsCDx85cr3V/vpJNOitr03Kj7Bn388ceJ7/fUU09Z7M+net/j9/Po0qWLxVnnW7nNS91jSPd2CyG+tujfDCHEe/GprCXTfb+0ktNZFXsPyfrIX6v2339/i/0x1/HVe6AXX3wx6pf0t70vya3lxbt16xa16Z5VukecXmc9P4a6H2dt7h/GShsAAAAAAIAc4qENAAAAAABADtWb9Ci/1Ol7folS2hK3pH6+jJcu79Yyjn6JpJacrkn6FUJo27atxZdccknUpsvYFi9ebPGvfvWrqF/WJWgsS/xBXf7/9bhreed27dpF/XSMfcnvch+vtdGmTRuLfYlvPa4zZ8602C8zXrNmjcU6njUpZ6npb5pSes4550T9dCm0xiGE0L9/f4u1bHEpfj/8/11LxWopyRDicpdZS5ZqymFNjp+WG9dUNk33CCGExx9/vKD3Rzq9B+rbt6/Fer0MIR7fUaNGRW2UYK+eLsu/++67o7Y+ffpY3LVr16hN5+nQoUMt9uc5PR+uXr3aYp9ir+k0ev0MIf4u6Jhef/31Ub+09OJyTSHX/9/GG28ctekxSUuJUnou9NdZ5f+G6dSpU5W/Fz8444wzLNaUwBBCqKiosFhLptdE0nfd/3va33eFzJdyKgeu/9eWLVtGbXouTDvmxx57rMXbbLNN1E+/F507d078XdqvWbNmiZ9RnwdoWrmn5+4Q4lT1rNs6FIKVNgAAAAAAADnEQxsAAAAAAIAcym16lF+WlLS828u6zFD7+eVWhx12WJX9LrrooqifVlNB9XRMtfKTr8Kgy+J0yZkuh6yJrEsPi7Fk0S+BLfW0uawpGQcddJDFupw4hBCmTZtmMVXYimfvvfe22J9P9RyqFTQ0HcrT9/DLwDVVMe38rOdMn+rj04KUX+payvz80HPKZ599FrWljZcqpIqTP5dpxQU9r/lqDn5cURw653TO+rmt34mHHnooaivlZfjFosfIV1/SNE2f1q2VUBo2bGixPx9q9btJkyZZrNfBEELYa6+9LL7ggguiNl3CrxX/fLWwNEkVcUr9O6Lnrtdffz1q03PcbbfdZrG/Nmnq1JAhQ6p87xDi9A9P06O0qs7777+f+DNeKY6V3jeeeuqpFuucCiFO//ZbJ9x6660W65jkJS2pFMctif5f582bF7XpvcOBBx4Ytel4a1rqkUceGfXLuv2Fnu+yPifw5+5FixZZfPLJJ0dtM2bMyPQ51hYrbQAAAAAAAHKIhzYAAAAAAAA5xEMbAAAAAACAHMrVnjZpuWmal5h1Hw3lf0bz2/r16xe1aY747NmzLX700UcTfy+qt99++1n8y1/+0mI/ZkuXLrX4jjvusNjnFxZSJjFrGfhClfoeNiFkP05asnTgwIEW+z0YJk6caHGh+1WVU45wEn98dL75ti+++MJin2ecRL/bX331VdSWNZd4hx12sLh169aJP+P3YNHvSKmPtd8rTffAaNy4cdSmpYDnzJljcdZrpN+3Rtu0fGYI8X4eOof9XhxZz4HM5+9k/b/qvhrbbrttYr8lS5ZUGaPm/PVI9xs57bTTorak8rD+e67zI22u6P4JgwcPjtr0HmnYsGGZ3i9tP49yuG+pyg033BC9HjRokMVaIvjee++N+unxWrVqlcVjx46N+ukY6j5UIYTQtm1bi7VU+3HHHRf1K2Q/svos7fqkdL7deOONUZvOl+nTp1u8YsWKqN/f/vY3i3VPEt0zKoQQttxyS4v9eOh1V++LyukalpXfl1T/LujZs2fUduKJJ1qsez61atUq6qfzT+eb7k8VQgg77bSTxWmlvHV/xptuuilqu+666yxeV3v3sdIGAAAAAAAgh3hoAwAAAAAAkEO5So8qpOygX/Kpr9OWfDZv3tziU045JWrTkplaavHrr7/O9Jmq+lzfK6clc7507YQJEyzWcm5+uaGWdNPlhmlj7SWV+Sun41/X/Hho6sb2229vsV9WqGX/so4P4/hjfimxlqAt9Dypsvbz817TtHT5+EYbbRT10zHVZa4hhDB37txMv7sU6HLfEEJ49dVXLT7kkEOithEjRlh8//33W7x8+fKon17vevTokfi7NZ3m7LPPjto0nU2/C++++27Uj7lZHH7O/vznP7e4WbNmFvvj/eSTT1qcVnIYa8cf97VNY/Hnb01H9Okad911l8U+nTKrck2JUr6ku5ZP1/QMX25aPf300xaff/75UZvOvy+//DJqu+aaayw+6KCDLN56662jfv78Wk6eeOIJi3/9619HbTpffLqLpo9us802ie9/5plnWuxTvtUGG2xgsT8vr1y50uJnnnnG4rPOOivqpynp5cqfM/Vv7SlTpkRtr7zySpU/l3Z/ofeUI0eOjNq6d++e+HN6DtXUuvHjx0f98pCqyEobAAAAAACAHOKhDQAAAAAAQA7x0AYAAAAAACCHcrWnTVZpJeGS8nR9TvBtt91mcbt27aK2t99+2+KXXnrJ4prk6pPX/+M9EVq0aFFlP83BDyEu0afHMa1kZdrxZizWTtpxT+u34447Wqw5wa+//nrUz5dhrE2lvL+RllgPId6DxI+Njkda+cOsNOdfc/VDCOE3v/lN4mdUmlN+8sknR21ahrHU+b3TdP+KTTfdNGrbZ599LL766qst9tc7LV2ssZYwDiG+FvrSmkr3apg5c2ZivzRpe5IVsudSfZH1fOr30Tj++OMt1jnrvy+33367xaV27OqLrOXstZ/uWRRCCBdeeKHF/vz33HPPWZx1jEvtepcm63Xel3TXa9W4ceMs7tevX9Rv8uTJFuu+K6tXr078Xffdd1/0WvfO0PPueeedF/XT++g87KlR2/T7fMUVV1i8ySabRP323Xdfi3WPrxDi+ww9V/q/F/V106ZNEz9T2neoTZs2Fg8YMMDi+fPnR/2uuuqqTO9XrgrdI0znuu6Heswxx0T9dKz9Xm96z6rz3p8f8oCVNgAAAAAAADnEQxsAAAAAAIAcqpfpUSpt+ZIumzrqqKOitmOPPdZiv7z4kksusbicluUXgy5L1OMYQjweuvRtyJAhUb88ln3WpXVt27aN2pYtW2ax/z6W0zJIv/S0d+/eFusS1TfeeCPql3UOF0rfQz+HH5s8LoWsCf//yfrd02XBvuylptnoEuQtt9wy6nfxxRdbfMQRR0RtSSlRflm/Lkt9+eWXq/nUpcsfFy137kvKdujQweL27dtb7Jf/6nXs008/tfjjjz+O+jVu3NjiLbbYImrbfffdLdZSnYWWHNZ5ueGGG0ZtaWkG9V3avNRjouVLQwhhl112qfJnFi9eHL1euHDhWny66ulnTEtDKYVrn7/+FHI9SjsOmgKnKfshxHPCl56dM2dOpvcv5XRgT/+velx9mkVaOpmW5f773/9u8bPPPpv4M1nT0yoqKqLXl156qcU33XSTxQceeGDUT0tWz549O2or9THVctonnHBCYj9N9w4hTu09/PDDLR40aFDUr0ePHhbrfZA/rprC788Bmral95c+pU5TcOr7vea65I9/r169LB41apTFPr1Y+ftLTSnO+9iw0gYAAAAAACCHeGgDAAAAAACQQ/UmPaqQZZ6bbbaZxSNHjozadNn/1KlTozZdOlXqyw+LTVOHmjdvnthPU9I6d+4ctekyUh0nPxa65N8vUU0aN78Mv1OnThb7aimaHrDrrrtarDuUhxDCo48+arH/Lk2YMKHaz5R3WT+338FfUxB1qeKMGTOK88EySkqP8t+ZvC+LrI5PidHUl6222irx53T5sB+bs846y+LDDjvMYl9pQV/7ykVKv0sffPBB1HbllVdaXN/Hopj0e7p8+fKoTV+nVXFKun76ua39HnvssahNlyHrPEob7zQ6xj4dqj5WPSp2Kopec0KIl+Hr8dHUihB+nO6tktJ7Cv28WSs4lopi/3/1mqnpjSHE5/MbbrghakuqrJKWzlXq46Pp2XpO8mnbOj/SzjPFHmv/u55++mmL+/fvb/EBBxwQ9dP57dN7NNW1FGUdA5/Wrfc+Y8eOtdhXX9RzrL6/T/l9/vnnLd5///0TP6P67LPPotfFSPVHCB07doxe6zzS9G5Px+P000+P2vz3J89YaQMAAAAAAJBDPLQBAAAAAADIIR7aAAAAAAAA5FBu97Tx+X9Zc0qbNGli8eOPP25xixYton5Lly61+Pjjj4/a0nLCkU5ztH3eteYZ6zhNnDgx6qe53Gk52fr+vjT7a6+9ZrHmNPfs2TPqpzmuvkRc0u/2+2307dvX4o8++ihqK/U8cnXwwQdHr1u3bm2xHrNCy1YWY98F/c7Ux30z0vjv5cMPP2zxBRdcELVpicw777zTYn+MtZ/Oo7Sx8G16nJcsWWLxOeecE/VbtGhR4nviO8U4n2QdOz8/dPy1jLu/tq5atarGn6nU5mKh9y/6c3369InadB8hvUY+8sgjBf2uQpXTNa0Y/9e0fWb22Wcfi33ZYt2DQffoKPR3lZq0/5veD6aVok97j2LsB5T2c7qP16233mqxv0ft0qWLxYceemjUNm7cOIvr074ctU2Pux6XF154Ieo3bNgwi/WapnEIIRx00EEW+70b9Xui95fPPfdc1K/UrnF1Sa99w4cPj9r8/cf3Pv/88+i17hs1f/78In66usVKGwAAAAAAgBzioQ0AAAAAAEAO5TY9KutyRF02FUIIZ599tsXdunWz2Kc8DR482OJly5YV8hFRBS3j+/bbb0dtWl5Pl9r78sHFWNKrpcd1eaRPgfJLkpUudVyzZo3Fb775ZtTv8ssvt3j69Ok1/7D1mM6/U045JWrTcdTSlAsXLsz8/sVOBynlUtL+WN1yyy0WH3XUUVFb165dLfZLgbNIS//wy7Tfeustiy+++GKLX3nllahfKY9NXSu09K/223nnnVPfE8Wnx/hnP/tZYpumoFVUVGR+/3JKbaoNafMo6/zQew5Nj9L08RDS75G0JLH+Xn8/XE4pGcW4fugx17jQ45g1FVXvKUePHh316927t8VbbLFF1LbhhhtanLWUeanIek3TY+Hvz2fMmGHxnnvuabH/uyDtHkm/d1OmTLFYS7X7z4GaadeuncU+bVi/BzoH9P43hBBeffVVi+vzWLDSBgAAAAAAIId4aAMAAAAAAJBDPLQBAAAAAADIodzuaZNVy5Yto9ennnqqxZqT6suvPfPMMxaT5108WsZQy+SFEJeE1tjnYWv+6Oabb25x9+7do36a560lUEMIYfLkyRaPGDHCYp933r59e4u//PLLqE3zWrWsqt8DScuXl1vZRc2p1vLpIcTH8+6777bYl2dH7dB9hI488siobeLEiRZvtdVWFvt9FZL2afD7B+i+GjrWIYTw5z//2WItY8t5t7jS9tQoZI+4Vq1aRW065lqevSZ7qpSyYnyfdf61bt068f11HrEX1LpRaFl37dejRw+L/X2Q3ttee+21Udtll11mse5vo/vwhVDa9yO1cf3QvS7S9r1IOtcW+pl0nMaMGZPY75NPPoleb7LJJhZrieP6vGdHEj8/9Fzpv/dJ50R/jz9w4ECLtXy67r/pf9cXX3wRteneNdddd13iZ0LN6N+BkyZNslj/3vJ0H9U//elPUVup/N3BShsAAAAAAIAc4qENAAAAAABADtXL9ChNWxk+fHjU1rlzZ4t1iUyyzTAAAAROSURBVOBdd90V9WPpWu3Q5aG6XDOEEMaPH19lXChdoqqpcCHU7ZJxTQkrN3rcfQnnmTNnWjxy5EiLSYupG3qc33///ahNUw379u1rcf/+/aN+ep7UdMR58+ZF/e69916L58+fH7VpGUbUnmLMK53PCxYsiNo+/PBDizXlrZRTMOpa48aNLU5LdZkzZ47F/tqH2lOMOab3JmnnRk3JOPDAA6M2TeV4/fXXLfZp4lnLIqNmavNY+nTTadOmWdyhQ4eoTb8jeh4oNHUvz3zKl772/9+k773vp+lmhxxyiMW9evWK+jVv3tzit956K2pbuHBhlb8La2e33XazuFu3bon9NF1twIABFq9Zs6Z2Ptg6xtUeAAAAAAAgh3hoAwAAAAAAkEMNarKcq0GDBrlY+6VViZ588smoTZcL6m7RPXv2jPq98847FteHJW2VlZXJpUFqIC9jWKamV1ZW9qy+W/XyOI4NGzaMXutu/zoX68N8S8NcLAklPRcLpXO4Y8eOUVvXrl0tfvnlly1euXJl1K8u53epzUVNcdhzzz2jtt69e1v84IMPWjx37tyoXz2sHFNWc1HHePvtt7f4zjvvjPq1aNHC4ksvvTRq02qommLl556+ru15WWpzMSn1KIQ4dbHY882n8GhqTpcuXaK2Zs2aWaxpOytWrCj015fVXCxV9XEu+upgWuFZ00P9eUwrofbr18/iengd9Kqci6y0AQAAAAAAyCEe2gAAAAAAAOQQD20AAAAAAAByqN6U/NY8zz322MNinwenNNf73XffjdrS8nspkwjUnC836l8DyDeds76su77mulg7NA9/6tSpUZvuI8Txr790jGfNmmXxvvvuG/XT+9AS2J+hXtB9bJo2bWqxlhUOoXbnn3/vVatWWez/jmnUqJHFfm8xoD5p0qRJ9Hrbbbe1WOfE8uXLo36nnnqqxeVwnmSlDQAAAAAAQA7x0AYAAAAAACCH6k16lNJl2lpKOIR4efeQIUMs/vbbbxPfz5fYY+kxAKCccR3MF8ajtKWV60bd0L8T8phutHr16sTXfF9Qn/nv9gMPPGBx3759LT7xxBOjfhUVFbX7wXKGlTYAAAAAAAA5xEMbAAAAAACAHOKhDQAAAAAAQA41qEkeZIMGDXKRNKl70Ky3XvzcSUt+lVKOZ2VlZYPqe1UvL2NYpqZXVlb2LMYbMY7rDnOxJDAXSwBzsSQwF0sAc7EkMBdLAHOxJFQ5F1lpAwAAAAAAkEM8tAEAAAAAAMihmpb8Xh5CWFgbH6QmNO0prZR3CelYxPfKxRiWKcax/mMMSwPjWP8xhqWBcaz/GMPSwDjWf4xhaahyHGu0pw0AAAAAAADqBulRAAAAAAAAOcRDGwAAAAAAgBzioQ0AAAAAAEAO8dAGAAAAAAAgh3hoAwAAAAAAkEM8tAEAAAAAAMghHtoAAAAAAADkEA9tAAAAAAAAcoiHNgAAAAAAADn0f7NB6mq1hghuAAAAAElFTkSuQmCC\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 19s\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.0015, PSNR (dB) = 28.21\n",
      "MSE = 0.0030, PSNR (dB) = 25.25\n",
      "MSE = 0.0009, PSNR (dB) = 30.23\n",
      "MSE = 0.0022, PSNR (dB) = 26.58\n",
      "MSE = 0.0035, PSNR (dB) = 24.61\n",
      "MSE = 0.0006, PSNR (dB) = 32.51\n",
      "MSE = 0.0051, PSNR (dB) = 22.95\n",
      "MSE = 0.0057, PSNR (dB) = 22.44\n",
      "MSE = 0.0029, PSNR (dB) = 25.33\n",
      "MSE = 0.0032, PSNR (dB) = 24.98\n",
      "Wall time: 11 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.9642\n",
      "0.9760\n",
      "0.9595\n",
      "0.9797\n",
      "0.9451\n",
      "0.9781\n",
      "0.9176\n",
      "0.9298\n",
      "0.9596\n",
      "0.9705\n",
      "Wall time: 3.94 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.9972\n",
      "0.9970\n",
      "0.9979\n",
      "0.9980\n",
      "0.9933\n",
      "0.9991\n",
      "0.9911\n",
      "0.9911\n",
      "0.9924\n",
      "0.9923\n",
      "Wall time: 13.3 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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Simple Autoencoders in Keras - 784-64-784-tf2.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
}