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

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import os\n",
    "#os.environ['KERAS_BACKEND'] = 'plaidml.keras.backend' # https://github.com/plaidml/plaidml/blob/master/docs/index.md\n",
    "\n",
    "from tensorflow.keras.models import Model\n",
    "from tensorflow.keras.layers import Input, Dense\n",
    "from tensorflow.keras.datasets import mnist\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# keras version still end up import tf.keras version which is incompatible, import directly instead\n",
    "from tensorflow.keras.callbacks import TensorBoard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoding_dim = 64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_img = Input(shape=(784,))\n",
    "encoded = Dense(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": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.compile(optimizer='adam', loss='binary_crossentropy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "60000\n",
      "10000\n",
      "x_train shape = (60000, 784)\n",
      "x_test shape = (10000, 784)\n"
     ]
    }
   ],
   "source": [
    "# Load mnist dataset\n",
    "(x_train, _),(x_test, _) = mnist.load_data()\n",
    "\n",
    "print(len(x_train))\n",
    "print(len(x_test))\n",
    "\n",
    "x_train = x_train.astype('float32')/255\n",
    "x_test = x_test.astype('float32')/255\n",
    "x_train = x_train.reshape(len(x_train), (np.prod(x_train.shape[1:])))\n",
    "x_test = x_test.reshape(len(x_test), (np.prod(x_test.shape[1:])))\n",
    "\n",
    "print(f'x_train shape = {x_train.shape}')\n",
    "print(f'x_test shape = {x_test.shape}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1595257665\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 45s - loss: 0.6934WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.169499). Check your callbacks.\n",
      "235/235 [==============================] - 15s 66ms/step - loss: 0.2433 - val_loss: 0.1583\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 18s 78ms/step - loss: 0.1412 - val_loss: 0.1253\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 20s 86ms/step - loss: 0.1174 - val_loss: 0.1080\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 18s 76ms/step - loss: 0.1037 - val_loss: 0.0973\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 17s 71ms/step - loss: 0.0948 - val_loss: 0.0901\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 17s 73ms/step - loss: 0.0885 - val_loss: 0.0851\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 16s 67ms/step - loss: 0.0842 - val_loss: 0.0815\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 18s 78ms/step - loss: 0.0812 - val_loss: 0.0792\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 19s 82ms/step - loss: 0.0792 - val_loss: 0.0776\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 17s 73ms/step - loss: 0.0778 - val_loss: 0.0765\n",
      "Epoch 11/50\n",
      "235/235 [==============================] - 16s 68ms/step - loss: 0.0767 - val_loss: 0.0757\n",
      "Epoch 12/50\n",
      "235/235 [==============================] - 15s 64ms/step - loss: 0.0760 - val_loss: 0.0750\n",
      "Epoch 13/50\n",
      "235/235 [==============================] - 16s 69ms/step - loss: 0.0754 - val_loss: 0.0745\n",
      "Epoch 14/50\n",
      "235/235 [==============================] - 39s 164ms/step - loss: 0.0750 - val_loss: 0.0742\n",
      "Epoch 15/50\n",
      "235/235 [==============================] - 32s 136ms/step - loss: 0.0747 - val_loss: 0.0739\n",
      "Epoch 16/50\n",
      "235/235 [==============================] - 26s 112ms/step - loss: 0.0744 - val_loss: 0.0736\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 25s 108ms/step - loss: 0.0742 - val_loss: 0.0734\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 26s 111ms/step - loss: 0.0740 - val_loss: 0.0733\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 27s 115ms/step - loss: 0.0739 - val_loss: 0.0732\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 29s 123ms/step - loss: 0.0737 - val_loss: 0.0731\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 29s 123ms/step - loss: 0.0736 - val_loss: 0.0730\n",
      "Epoch 22/50\n",
      "235/235 [==============================] - 25s 105ms/step - loss: 0.0735 - val_loss: 0.0729\n",
      "Epoch 23/50\n",
      "235/235 [==============================] - 31s 133ms/step - loss: 0.0735 - val_loss: 0.0728\n",
      "Epoch 24/50\n",
      "235/235 [==============================] - 31s 131ms/step - loss: 0.0734 - val_loss: 0.0728\n",
      "Epoch 25/50\n",
      "235/235 [==============================] - 35s 148ms/step - loss: 0.0728 - val_loss: 0.0723\n",
      "1595257665\n",
      "Wall time: 24min 45s\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoded_imgs = encoder.predict(x_test)\n",
    "decoded_imgs = decoder.predict(encoded_imgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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+WH4yJAAAUQ0lEQVR4nO3beXCN59/H8Ss5EYQTiVhjiX0pYg8ZP0vtzKBqRnVqqRpFO8wYY1+qythqGEtrq2EGU0trm9r3XSy1ZchIrUUkoWQX4Tx/PH8+nF6fmeP5Xe28X3+ad7/u5Jz7Pvf9dRrk8/kMAAAAAAAA3BL83z4AAAAAAAAA/F8sbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHBSixKGhob6iRYta9wUFBdLBlC1bVurz8vKs22LFikmzHz58KPXR0dFSf+fOHev2zZs3xufzBUl/wTsUL17cFxUVZd2npKRI8z0ej9SXKlXKun39+rU0Oy0tTerr1asn9VlZWVKfnJyc7vP5Skv/0TsUKlTIV7hwYeu+XLly0nxltjHa7zokRLrsmCJFiki9+h4MCrI/tZ48eWJevHgRkHMxNDTUp/xs6s+lvj+Va1hYWJg0+/nz51IfHKz9e8KrV6+kPi0tLWDnYmRkpK9ChQrWvXLtN8aY0qW1w1Tez+rrmJqaKvWRkZFS/+zZM+s2KyvL5OXlBeRcLFy4sE/5XSifW8b872e4Ijc317r1+XzSbPXcCg8Pl/qMjAypf/ToUcDOxRIlSviUzzr1WNXXMScnR+oV6vW9atWqUq+ci7m5uebly5cBORe9Xq90j6pek9T7D/XeSXHv3j2pV++FKleuLPXXr18P2LkYERHhU+4plGueMcakp6dLfcWKFd/bseTn50t9dna21CvX7JycnICdi+rnovqMpn52FSpUyLpV7z0yMzOlXn1GUt+veXl5bz0XpatX0aJFTatWraz7p0+fKuPNyJEjpT4pKcm6jYuLk2ZPmzbtvfaDBg2ybtUT3J+oqCgzadIk63727NnSfPVEGTJkiHX7119/SbNXrFgh9UePHpX648ePS33Pnj21T2g/ChcubBo0aGDdjxs3TppfrVo1qV+1apV1q75H6tatK/XqQ4ZyI/TVV19Js//u71WuSxEREdL8kydPSv2ECROs26ZNm0qzt23bJvXve8m+fPnygJ2LFSpUMFu2bLHuBw4cKM1X33PKQ0njxo2l2YsXL5b6Pn36SP2mTZus2507d0qz/QkLCzPt2rWz7ocOHSrNVx/er1y5Yt2q/zjm9XqlvlOnTlJ/4MABqZ88eXLAzsVy5cqZ5cuXW/eHDh2S5qs39pcuXbJu1WWCer8yY8YMqd+6det7OxZ/oqKizDfffGPdq9ck9f5j4sSJ1q360Dp8+HCpr1OnjtQvXbpU6mvWrBmwczE6OtqsX7/euk9MTJTmr1y5UuoXLFhg3V6/fl2arS7fzp8/L/XKPeqRI0ek2f6EhYWZDh06WPfKotcYfQmu/ANW3759pdmHDx+W+urVq0v92rVrpT4xMfGtbyr+9ygAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB4VIcUiIKVmypHW/e/du6WAOHDgg9R07drRuR40aJc1u1qyZ1J84cULqBw4caN1u3rxZmu1PQUGBSU1Nte6vXLkizR8xYoTUly1b1rrNzs6WZvfs2VPqFy1aJPWRkZFSH0her9d06NDBuh87dqw0v1q1alJfq1Yt67ZEiRLS7Pz8fKn/7bffpH706NHWbWhoqDTbn4iICNOjRw/rvkuXLtL8DRs2SL3yfrp69ao0u2rVqlJ/8uRJqS9evLjUB9KdO3fM559/bt1nZmZK88PCwt5bP2/ePGl2TEyM1Ofm5kr9xo0brVv1uuBPsWLFTHx8vHW/b9++gP3db3Pr1i3rVrmGGGNMcLD2b3Xbtm2T+vHjx0v95MmTpd6f4OBg6VpQpEgRaX7btm2lvmHDhtbtzz//LM3evn271FeoUEHq7927Z90G8lwsVaqUGTx4sHV/6dIlaX5sbKzUL1iwwLpV7w+8Xq/UJyUlSf3cuXOlPpCeP39uduzYYd2rzznNmzeX+rVr11q36rNr165dpX7YsGFSr5xf6jObPx6PxxQrVsy6T0hIkObXr19f6mvUqGHdKuetMfozv3q97ty5s9QnJia+9c/5pg0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINClPjZs2dm06ZN1v3HH38sHcz8+fOlvnv37tZtVlaWNDshIUHq69atK/W3bt2ybjMyMqTZ/hQUFJj09HTr/tSpU9L8OnXqSP2JEyes2+TkZGn28OHDpX7t2rVSX7p0aakPpMzMTHPkyBHrftSoUdL8nTt3Sv3Lly+t27Fjx0qzmzRpIvUxMTFSP3DgQOv27t270mx/ihYtamJjY6179dwqWbKk1E+aNMm6TUlJkWZPmTJF6g8fPiz1zZo1k/pACgsLMw0aNLDuExMTpflr1qyR+qioKOs2Li5Omt2nTx+pnzlzptQvWbLEup01a5Y025+UlBQzZ84c6169t7l9+7bUd+zY0bqdPXu2NPvBgwdS/5///Efqz549K/WBlJ+fb+7fv2/d9+jRQ5q/cOFCqW/YsKF1GxERIc0+f/681J87d07qGzdubN2q76m/m6Xcr5w+fVqa7/V6pf7ixYvW7aBBg6TZxYsXl/rNmzdL/ZYtW6Q+kAoKCszTp0+t+969e0vzlfenMdozoPKaG2PMzZs3pV69N1Pe46mpqdJsf/Ly8qTnLvU+7NixY1KvPDd369ZNmq1eR5o3by71e/bskfp34Zs2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOClFir9drWrRoYd1nZmZKBxMfHy/1GRkZ1u3cuXOl2b/++qvUR0dHS/2xY8es2z179kiz/QkPDzddu3a17hctWiTNb9mypdQfP37cuh08eLA0e/Xq1VLv8XikPiEhQeoDKTs725w5c8a6T0lJkeb37dtX6nft2mXdrly5Upo9bdo0qc/Pz5f6mjVrWrdJSUnSbH8eP35sZs2aZd0nJiZK89esWSP1n376qXUbGxsrzU5PT5d69VzfvXu31AeS1+s17dq1s+7V69K2bdukfsqUKdbt77//Ls3OycmR+i5dukh97969rdsVK1ZIs/3xeDwmPDzcur948aI0v2nTplK/bt0663bQoEHS7IULF0p9SIh0m2gePXok9YH06NEj6f0/f/58aX5ERITUP3z40LodMGCANPv8+fNSr35+KOfuoUOHpNn+BAUFmcKFC1v37du3l+bfuHFD6pX7iWrVqkmz9+7dK/VDhgyR+tq1a0v9mDFjpP7vBAUFWbfdu3eXZt+6dUvqT506Zd0WFBRIsytWrCj16mf6jBkzrNsJEyZIs/3x+XwmLy/Pulffb5cuXVIPydrixYulfvTo0VKvPLsaY8zNmzel/l34pg0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOChEiV+9emVSUlKs+zlz5kgHM3LkSKm/dOmSdRsbGyvNjouLk/rs7Gypb926tXXr8Xik2f5kZ2ebs2fPWvf79u2T5sfHx0t9ly5drNsHDx5Isxs2bCj1LVq0kPpVq1ZJ/eXLl6Xen9jYWOm1KVeunDR/ypQpUv/hhx9at1FRUdLsL7/8UuofP34s9Xl5edZtcHDg9txFixY1DRo0sO6V89YY/f3cv39/63bhwoXS7L1790r9jh07pH7mzJlSH0h//vmnGT9+vHW/ZMkSaX6vXr2kvlatWtZtcnKyNDssLEzq+/TpI/WtWrWybm/evCnN9icmJsYsW7bMuu/UqZM0/4svvpB65femfIYaY0xCQoLUe71eqS9fvrzUJyUlSb0/oaGhJiYmxrrfv3+/NL9ly5ZSv3XrVuv27t270mz13Fq+fLnUK79H5TP076Snp5s1a9ZY97m5udJ85ecyRjt3jx07Js3OyMiQ+kqVKkm9z+eT+kAKDQ2Vjnf9+vXSfPWZTumrVKkiza5du7bUf//991KvXlMD5c2bN9K5rd7nqe/nyMhI67Z9+/bSbPW54dmzZ1KvXq9/+eWXt/4537QBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHBQiBJHR0ebqVOnWvdr166VDubp06dSn5WVZd22bdtWmv3gwQOpr1+/vtSHh4dbt8HBgdutZWVlmdOnT1v306dPl+YPHDhQ6lNSUqzbvXv3SrNfvnwp9ZGRkVJftmxZqQ+kBw8emFGjRln3Z86ckeY3adJE6pVj6dGjhzS7atWqUh8fHy/1jx49sm5fvXolzfYnIyPDHDx40LqvWLGiNL9169ZSv3TpUut22LBh0uwLFy5I/cKFC6V+9OjR73W+PwUFBSYtLc26V6+R+/fvl/pGjRpZtzVq1JBmHz16VOrHjRsn9f369bNuP/vsM2m2P8+ePTMbN2607gcPHizNj42NlfouXbpYt0OHDpVmN27cWOojIiKk/vnz51IfSDVr1pTuE+bPny/N37Vrl9SHhNjfYt+7d0+arfbKuWWMMQ0bNrRut2/fLs32x+PxSPfHPp9Pmt+3b1+pT01NtW7v378vzR4yZIjUr1u3TuorVaok9YGUk5Mjfe5XqVJFmq8+Dyj3bnv27JFmf/vtt1J/5coVqT9y5Ih1q96b+RMeHi59Ft26dUua//r1a6lXrqc5OTnSbPVzrkOHDlLv8Xik/l34pg0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINY2gAAAAAAADiIpQ0AAAAAAICDWNoAAAAAAAA4iKUNAAAAAACAg1jaAAAAAAAAOIilDQAAAAAAgINClDgzM9McPXrUuo+Li5MOZv369VJfvXp16zY5OVmaPXr0aKn3eDxSn5qaat0GBQVJs/158+aNyc3Nte6fPn0qzW/QoIHUP3z40Lpt1aqVNDs6Olrqx4wZI/WdO3eW+kAqUqSIqV27tnX/+PFjaX7Hjh2l/uLFi9ZtmzZtpNlt27aVevW6M2fOHKkPlDJlypivv/7aur9//740f968eVLfqVMn63bmzJnS7JiYGKlXrlHGGJOUlCT1gfTBBx+YzZs3W/dLliyR5r9+/VrqleuS8nlujH7uvnjxQurPnTtn3WZnZ0uz/QkJCTFRUVHWfWxsrDS/efPmUr93717rtqCgQJpdo0YNqR86dKjUz549W+oDKSsry5w8edK6P336tDR/yJAhUp+enm7dqp+506dPl/rz589LfZ06daxbn88nzfYnNDTUVK5c2bq/ceOGNP/atWtSv23bNuu2f//+0uz9+/dL/ciRI6X+jz/+kPojR45IvT/h4eGmW7du1v327dul+WXKlJH6iRMnWrfq66jeo96+fVvqV69ebd0q15y/8+bNG5OTk2PdZ2ZmSvODg7XvjSjnrvq8qB5LSIi0PjFXr16V+nfhmzYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA5iaQMAAAAAAOAgljYAAAAAAAAOYmkDAAAAAADgIJY2AAAAAAAADmJpAwAAAAAA4CCWNgAAAAAAAA4KUeKSJUuafv36Wfc//vijdDAfffSR1KelpVm3t2/flmYnJydLfUJCgtSXL1/eus3NzZVm++Pz+czLly+te+V3bIwxly9flvrnz59bt8uWLZNmt2nTRuq/++47qZ86darUB1Lx4sWln69nz57S/PHjx0t9XFycdatcQ4wxZsOGDVLfokULqW/SpIl1e+3aNWm2P6mpqWbJkiXWfb169aT5n3zyidSvWrXKuh0+fLg0OzQ0VOqvX78u9SVKlJD6QLp586Zp3bq1dX/s2DFpfv369aV+xIgR1u3cuXOl2XXr1pX6H374Qep/+ukn6/bFixfSbH+8Xq9p3769dd+hQwdp/syZM6VeuSaNGjVKmq1ewx4/fiz1FStWlPpAys3NlX6+0qVLS/OV67UxxnTu3Nm63b9/vzQ7KytL6rt16yb1Bw8etG4zMjKk2f5ERESYXr16WfdKa4wxkyZNkvpGjRpZtxcuXJBmlylTRupXrFgh9QMGDJD6QMrIyDD79u2z7qOjo6X5T548kfpdu3ZZt8HB2vcZRo4cKfXK78UYY6ZMmWLd5ufnS7P9KVasmGnWrJl1v3v3bml+fHy81CuvealSpaTZHo9H6q9evSr16rP0u/BNGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAAB7G0AQAAAAAAcBBLGwAAAAAAAAextAEAAAAAAHAQSxsAAAAAAAAHsbQBAAAAAABwEEsbAAAAAAAABwX5fD77OCgozRhz7/0dDt4hxufzlQ7EIF7D/ypex38+XsN/B17Hfz5ew38HXsd/Pl7Dfwdex38+XsN/h7e+jtLSBgAAAAAAAP8/+N+jAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABzE0gYAAAAAAMBBLG0AAAAAAAAcxNIGAAAAAADAQSxtAAAAAAAAHMTSBgAAAAAAwEEsbQAAAAAAABz0P6xAcbvYEzxgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZQU1dnH8TtvXMENkE1FFmXABaIociSKC1ExURPjEgX3LS5oYlw4EhMxHqOJJigk7gouqHgkKhI1gom7BlQUg+LGoiLIqqCgxqTfP3Ly+LsPU0VN0zPUVH8/fz3DvVPT1O1bVd3nPvepKZVKAQAAAAAAAPnyf2v7BQAAAAAAAGBVfGkDAAAAAACQQ3xpAwAAAAAAkEN8aQMAAAAAAJBDfGkDAAAAAACQQ3xpAwAAAAAAkEPr1KdzTU0N9cHXklKpVFOJ4zCGa9WiUqnUuhIHYhzXHuZiITAXC4C5WAjMxQJgLhYCc7EAmIuFUOdcZKUN0HjmrO0XACCEwFwE8oK5COQDcxHIhzrnIl/aAAAAAAAA5BBf2gAAAAAAAOQQX9oAAAAAAADkEF/aAAAAAAAA5FC9qkcBAAAAAFBNamriwkylEgWW0HhYaQMAAAAAAJBDfGkDAAAAAACQQ6RHodH83/998x1hu3btLB4+fHjU77vf/a7FG220kcVfffVV1G/hwoUWH3HEEVHbK6+8YjHLF5sOv/RU3zP/+c9/LGZM80XHjbEBAABFt84633yM/vrrr9fiK0E1YKUNAAAAAABADvGlDQAAAAAAQA7xpQ0AAAAAAEAOsacN1ojfg2Tddde1uEePHlHb6NGjLe7evbvF3/rWtzL9LT12CPF+N88991zU1q9fP4snT56c6fhYOzQnWPczCiGEM844w+Jx48ZZPGbMmKjfv//97wZ6ddVN5+Yvf/lLi88999yon+49dPDBB0dtTz31lMVZ97vR4/nfY8+cb+j1d7311kvs969//cti3RsKxebvz0l7TzGnKkuvm/65RffmYy4CjacS5br9s0nLli3rbFu6dGnUT+/BzZo1SzyGvka/j+e8efPq/XpRLKy0AQAAAAAAyCG+tAEAAAAAAMgh0qNQUT179rT4gQceiNq22GILi3UZoS+TN2XKFIuvuuoqi/v06RP1O/vss+s8XgghtGnTpj4vGzmx8847Rz937tzZ4uXLl1vMcv6G4efRQQcdZPGFF15o8YYbbhj106W/tbW1UZumR6XRZcGaMhdCnP5Wzalwfnn3+uuvb7GmnDZv3jzq989//tPiZcuWWVzuPNL3SVp6qx8r0kGyK3cpv47NAQccELUNGTLE4nvuucfi2267Leqn8xnZbLnllhafcMIJFn/wwQdRvz//+c8Wf/bZZ4nH0/H312WdR2nvCz1GixYtorb27dtbPHPmTIu/+OKLqB/32lXpedXrn7++Zb3eJaUt1uf3NOY6m03SHNP7aghx+tKBBx4YtZ155pkWd+jQwWKfFrl48eI6/24I8ecVfR0+PUo/A82YMSOg+rDSBgAAAAAAIIf40gYAAAAAACCHcpse5ZePKZZr5kfaOPkl1rpUXpfoDxw4MOr32GOPWaxj/cILL0T9dFmiX4p4+OGHW/zoo4/W+RpQWeXOWV0O6tNudAf+qVOnWszy34bRrl276OfLLrvMYl0y7Mfzyy+/tPidd96J2nR80+afHtMvC8Z/+fOuc04r5rVq1Srq9/bbb2c6ftoy/aSl5D4VVefwhx9+GLX51AtU3iabbGLxTTfdFLVpCk/Xrl0tvvvuu6N+pEetXpcuXaKf9Tljgw02sPiKK66I+um1Mqusz7z+Hty2bVuLL7nkkqitU6dOFl9wwQUWT58+vd6vr4j0GqdzJYT4mVXTVPyWAGnXu6TrqR/rtGcdPgutGf3ccN5551k8dOjQqJ/OZ5+q6H9OotVu/ZgmpRh//vnn0c++6hSyS5tjTWkesdIGAAAAAAAgh/jSBgAAAAAAIIf40gYAAAAAACCHGn1Pm/XWW8/iXXfdNWo77LDDLN5uu+2ito8++sjicePGWfzWW29F/XQPDM0Z9XskaJlpnwecVEbP0xy5FStWRG3VunfKnDlzLB4zZkzUpuXqtJT3E088EfVLyi/s379/9LOWtfW/85e//GW1x8Pa4eeUlkk85JBDojbNr583b15Zx1e8F/5Lr11aVl3nTQjxvg1p51X3MRk7dmzUpnN91KhRFi9ZsiTqxz5F9dezZ0+LzzrrLIu1xHcIcWnhcvfH0N/TsdLSwSGEcOyxx1o8fPjwqE3vD8zFyvDjNGDAAIu32GKLxL7vv/++xf75BXXbeOONLda990IIoWPHjha/9tprFmuJ7xDiZ8805dzHfKnim2++2eIddtghapswYYLFui9LNc9LvS/q55Hbbrst6qf3u4ULF1rs9wPS67C/v+l5rsQ5r+Zx85L2mfH/rtfHiy66yGLdf8ZL23NIY39N1f3dXn/99ahN98zRz7CXX3551O/dd99NfF3VIu3zul6f99tvv6jfiSeeaLF/ZtHvB6ZMmWLxjTfeGPWbPXu2xX6/Kn1fNOTnf1baAAAAAAAA5BBf2gAAAAAAAORQTX2W1NXU1JS1/k6XL7Vs2dLiESNGRP32339/izX1JYTkctHLly+P+mm5S+X7aak3vxROl5KnlZ7VlC0tmRhC5csTl0ql5LWy9VDuGNbj+Bbrkr8Q4nJ1n376qcVpy4U1ne69996L2rR8qU+T22mnnSxOG8NGXlL6cqlU2nX13VavocexHFmXc/vyhueee26dcQghnHLKKRYnlYJf3euo9Bg3lbno/lb0sy7l15So2traqJ+OVdr4Kn++tXzwypUrLR40aFDU75FHHkk8RgNoknNxnXXijGYtMdu3b1+Lzz///Kjf7bffbrHej9LKl/olvjom+jp+85vfRP2OPvpoiy+99NKo7dZbb63zeOVqinMxq6zXU71HhhDCxIkTLd5zzz2jNr3X7rHHHhZPnjy57NdZAbmdi/qcGEIIV199tcUnnXRS1LZ48WKLv/Od71g8d+7cxOPrfEu7b6WVqNVjHHjggVG/O++80+KPP/44attnn30snj9/fuJrzKoIc1E/P7z88ssWb7vttom/o/c3LfseQvw8o6kVIeQ2HTi3czHNZpttFv2szzf62cB/FtDPo+PHj7e4R48eUT/9zKnbc4QQwuOPP27xJ598YvGbb74Z9dPPPFnHvtz3SBHmotLPkn4uDhs2zGJNifKfP9OutUnpdL7k+qJFiyx+5ZVXojZNBX/++ectXoN5XudcZKUNAAAAAABADvGlDQAAAAAAQA7xpQ0AAAAAAEAONXrJb93TwJfT0pJobdq0idp035lWrVpZ3KJFi6if5qTq39L9VEJILxeteWyaL96pU6eon5YqPu6446K2adOm1Xn8opfl0/+fnv8QVi2RloXm5Ldr1y5q0z0XfvWrX0VtmmesrynrvhzIppyypH6fAM3D13ELIXkepb2Oht7Tpiny+whdd911Fnfv3t1in9ubVObZn1Od634M9dqtJRl1v4UQQth+++0trsQeC0WkJdhDiPfu0n3bHnrooahf1rzqrPcqncO6F10I8b1Vc8BXd0xkv55qP/+stPPOOyceT8sTa1nq+rymoo+h/n87d+4ctem9Sve5CCGEY445xuK0fWySpJWETnuNm2++ucV+fyk95jXXXBO1+T1uqpF/b/fv399i/7yf5RgDBgyI2nr16mWx/4zw1FNPWZzT/W1yTcfn4Ycfjtq09Lrui+g/gyxZssRiHXst6R5C/JnQP98kzdOiXycrzT+jduvWzeK7777b4q5du0b9dKx0HunetCGE8MQTT1js75n6HOVfh9IS8b5suP6e3gv861hTrLQBAAAAAADIIb60AQAAAAAAyKFGSY9KSpl57rnnon5aJssvLUsqJb3ppptG/XTZti6F0xJwIcRLr5588smoTZe96pLG3/72t1E/TSV49913o7a0VIJqlfU8aCqbLun1y9a+/PJLi33J0qSl5IxF3RpzCbxPc+vdu7fFPp3CL0HPgqXGq+rZs2f0sy4FTip3GEL8PtCx8OUOf/3rX1s8b968qO33v/+9xZpKo6lSIcTpsgMHDozafOnFarL++utbPHjw4KhNx+eSSy6xeOnSpZmOnVZKOI0uL27btm3UpmVPn3322UzHK5pKXE+zpsccccQRUZumI3o6Hn6ZP/5Lz22fPn2iNi0t7K9zU6dOrfffKvdepSkBF110kcX+3vr0009bfNddd0VtWd+T+tylKelF4OeKXkP1faDPmiGEMGfOHIt1C4Wtt9466rfVVltZrCkeIYSwxx57WDxz5kyL0z771KcsfNH4sXrssccs1rSVEEL461//arEv8610/ulnU7+1A58hGoaOqc69EEI46aSTLNaUa3/N1OeNkSNHWnzVVVdF/XRM9ZkqhBD69u1rsX5XcOqpp0b9NN3Rb/PQr18/i/X5yD+7run7h5U2AAAAAAAAOcSXNgAAAAAAADnU6NWjdGlQuUstdTfmtKVHuqRtwYIFUb8pU6ZY7Jdb6e+tWLHCYp9GoMuy7r///qitaMtIK03PsU+TuPnmmy3WijJ+WdmkSZMs9kuVWcLY+NLOuc4d3Vk9hHip95tvvhm1JVUcowrY6q2zzjeXd60WFUK8pFv5a+GsWbMsHjVqlMXXXntt1C8tfUlTeh5//HGLt9lmm6jfAQccYLGmW4UQpwCkLXcuAv/e7tGjh8W+UtM777xjsd6Dss7F+qRHaZrEOeecY7GmjIQQVyXSeyRWL+t9S5d36zLyEOL3z9dffx21aapiORXFqoG+zzWNNIR4mb5PydBKU2+88YbF5Z4/HUdfJXXYsGEWDxo0yGKfxnPxxRdbrM+y9VG0Z1k9r4ceemjUpvckvXZp+kQI8XYItbW1Fr/66qtRP52nrVu3jtqOPvpoi6+44gqL/fmutvmn9F41bty4qE3Tcn269uWXX26xpoGWm7aa9XkzrVJtNY/j//iqXPpZ7/DDD4/adOz1OeeMM86I+r3wwgsW6/Uv7Xz7Z0idt/p+8duq+JQopXPdf6atJFbaAAAAAAAA5BBf2gAAAAAAAOQQX9oAAAAAAADkUKPvaVNpaXlr5ZbD09y0888/32KfA657PMyfPz/z8RGf4zvuuCNqO+iggyzWc64l/kKIc4KLvs9FQ2jMMoa6h8rBBx8ctem+NX/84x+jNr8nw/+QH7x6Wk5xhx12SOyn5/KZZ56J2k4//XSLtSxpfebb3LlzLb711lstvuyyy6J++h45+eSTozYt31j0a63Pmz7ttNMsbtmyZdR2zTXXWJy2Z4Xmhzdr1izxd9Jy8rt3726x7qPi93p78MEHLS7afhhZNcT1ScdD93Po1KlT4u/onA0hhNdff73ir6vIli9fHv2sY+D3cpowYYLF48ePt9iXAl+0aJHFeu/zJY11L6vjjjsuatMx1/faQw89FPV7++236+xXzXTPohNPPDFq02uv7qum59HTNv9+0edcfz3t0KFDxlf8jWobw912283i3r17R226t6k+p4QQwtKlS+s8nr9X6ZjoZw0/VlmfldP2i6NseHzfCiHex1D3YAwhfnbQPRSfffbZxH56jv1Y6/W6W7duUdvuu+9u8ZFHHmmx34dK+TGcPXu2xbo3Z6XHmpU2AAAAAAAAOcSXNgAAAAAAADnUZNKjkkquVWLpkS6XDCGEESNGWKwlABcuXBj10zJ91brcrT5l7bSvnldfxlaP8dJLL1nslwinlRnWv5W2LLFaxy2Eyvzfsx5Dlxn6Us9arv35558v6/hYVbt27SzeYIMNojY9r48++qjFP/rRj6J+WkKx3DLrmuKm6VEXXnhh1K9Vq1YWa2ndEEL49re/bXHR06N8CpQu3fUeeeQRi7POFU3JyFr2OYQ4bVVTOXSZeghxigbzt2FoKWpfRlWXi99www1Rm4496qbXq7Fjx0ZtupxfS3yHEMKWW25psaZr+Dmgzy1aXtanKmpqjU8r0KX/K1eutPjSSy+N+iWlF1czvY+1b98+atProS8jncR/fshK0zXKvbcWMeVG06T/8Ic/WLzppptG/XTu+M9mSbLe7/y5zJrmW4kS4kUZx7roM14Iq6ZEKT0Pffv2tdin5mvqf58+fSzWNPAQ4uuzT0VV+v5Lm9tLliyJfj7vvPMsbsjrLittAAAAAAAAcogvbQAAAAAAAHIot+lRfudnlVbhIiv9PV9ZRasSqeHDh0c/L1u2rKy/XST1WcqnS81Gjx5tsS4DDiHegX/o0KEW+6XdaWOflB7lFXkpYkMpZyf9H/7whxZvsskmUb+JEydanJbylvQaVvc6qpUuJ/bnS9Oezj777Dr/PY0uIQ0hXnactjRUxzdtybG//vtrRJH5/6suIf7kk0+iNq1Ek0bHJ+sScf+e0fukjt2LL74Y9cu6VD3tb1XTfM56PdU5d/zxx1vsl3Brutq4ceOitmo6r+XSc+Tf2/vtt5/FgwcPjtp0ab6mpvpz/tFHH1k8efJki32lL62cM2jQoKhNrwk6xu+//37IKukZqT4pk02R3p/02SOEuGrhPffcY/Hf/va3qN+MGTMs1tRdf19Mq7Ko75cuXbpYnFapqhpoKrdWSfPPBC1atLD4jDPOiNp064pyU6KSVPO9ak2988470c+aSn3EEUdEbTpu+vnhqKOOivrp/S/tOUf7pX2/kPQaQghhzpw5Fg8cODBqmzJlSuLvVRIrbQAAAAAAAHKIL20AAAAAAAByiC9tAAAAAAAAcihXe9pk3Xska1npNJrflpYjp7mrf/rTnxJfE1Zvr732sninnXay2I/Z1KlTLZ4+fbrFPndff/b7aCSNDWPWeHRvDi3X7sdg/PjxFqftc1Lu/lXVqnfv3hb7c677Q3388ceZjpf1/Kf123HHHS32JTyVlvMMIb4OFJ0v/av72Phc7M0339ziDz74INPxdXz88fTn2traqG3AgAEW6/tpwoQJUb+s5VHT9h3jOr0q3cNB99HwdE+hBQsWNOhrKjr/XKF7Glx44YVRm84d3XPGv7f12qZ7H/i5+NJLL1m89957R226N9hPf/rTOo+3OkUsF52F/l9/97vfRW26d5BeW3/84x8nHkP3kPJlwufOnWtxv379orY2bdpYfNNNN1ms+yaFsOq9MOl1FIXuA7R48WKLdZ+oEEJYd911LR42bFjUdswxx1j85JNPWvzhhx9G/R5//HGL9TnIP5vo3/L7/unv6etNm4tFHLcs/D6wp5xyisXXXntt1HbYYYdZrJ8Xu3XrFvXTOTZ79myLd99996hf586dM73GlStXWnz99ddHbfo+03kfQuONKSttAAAAAAAAcogvbQAAAAAAAHJoraZHpS2jL7csqUpbrqRl5b7//e9Hbboc/YQTTrBYl03VR7UuQ9VzHEIId9xxh8Wa2uTLPGvZxaVLl1rs3xNZl+Grajr/lVJu6kLLli0t3nrrrS2eNWtW1E/TK9KOzdil8+OkS719amE51ySdf2mlTHUpcQhxKsf9999vsS+PqnxpSL+suciWL18e/fzss89a3L9//6jtyiuvtPiZZ56x2JcN//TTTy3WZfnt27eP+mk6iC/B2bp16zr7+VLF5ZROraa5nfV66vtpumPz5s0t9vfBG2+80eK0eYo148dNx6ESzyY9evSw2KdOjRgxwmJ9Xk2TVua2muaf0vLrIYRw1VVXWTx06FCLmzVrFvXTFNaRI0darKWmQ4jTkP3njLvuustivUdus802UT/doqEa6DXrvvvus/jiiy+O+unzg3++0RSarl27Wuw/QwwZMsRinQNp5aI19TGE+N566aWXWjx69Oion74XqnW+eZpqpumgIcTbZOi90I+Nnkt97tG5HEJcFt6ff03b0lTISZMmRf3Kua5XGittAAAAAAAAcogvbQAAAAAAAHKIL20AAAAAAAByaK3uaVNuXp/mt/lcRl+i8X98Pq/mHnbs2DFq0zJwr732msX1eb2UJ473AwohLtmn51Lzs0OIy7Zpv3LPKfmjjcOPj5YI3mijjSx+/vnno36+xDHK4/eS2X777S32Y6PXzUrMDz2ejnsIce7+xhtvnPiaNJddS3aGsGqZzSLz/9cbbrjBYj/G++yzj8VaFtjn3es1VXPC/fF0vymds56+Z956663Efmn0npyHXPHGknW++TE855xzLNZzp3slhBDCvffea3F9SkCjcsrZa9GXhL7uuuss9uVlJ06cuNrj+dfh+xXtuWhN92kLIX4W/fvf/26xLzOsnwveffddi/3+jPo69HNFCPFeNd27d7fYl5I/7bTTLE76fFMkes6GDx9usb9X/eQnP7F4k002idqS3gv+86LuVZR1zyf/3tK93q6++mqLN9tss6iftqWVccd/lfNM0KVLF4t13oQQj73/zHHWWWdZrPvY5PG5hJU2AAAAAAAAOcSXNgAAAAAAADlUU8+Un1ysp0xb2q/LHXWJ3OGHHx71u/POOy32y4u/973vWfzCCy8k/q3GVCqVKpJv1dBjqEvqfdleTY/Sc77zzjtH/XS5vZ5zv3wxbcliVmnLmPXv6VJ1/ztZUzdKpdLLpVJp13q+xDrlZS4qv3z1ueees7hXr14W+3S4888/32K/XLkSZYErXVo4r3PRp7MsWrTIYl8CesmSJRZrOXa/bFSvtbqUWH8nhDh149hjj43aNthggzpfrx/r8ePHW3zkkUdGbQ2wnLjJzEW9Dvkxrq2ttVjLxvp0ijlz5lis194FCxZE/TRFbdSoUVGb3hf1fbLVVltF/bQEahq/VF1lXZac17lYCa1atYp+1jQMXXrv77NaKrqJpBU2ybmYNS3JS0pPnD59etRv2223tdin1mj56LS5Uul7X5q1PRf1/5pWIjjrMXSs/XgmpXamjYW/3p1++ukWazlrLT8cQnzdnTlzZtTWAGOaq7mYNKYhpD+T632yZ8+eFmu6VQgh7Ljjjhanfa7UtDc/jknPN9OmTYt+3m233SzW+2xDWNtzsSH587/vvvtarM+Qflx0TB944IGobeDAgRbn6J5Z51xkpQ0AAAAAAEAO8aUNAAAAAABADq3V6lHlyrp0ukOHDhbfeuutUZsurXvqqaeitldeecXiou2w39B0SW/Lli0z/Y7unB9CCO+//77FuuxR0zN828qVK6M2XX6oy2P9MfQ90qdPn6hNl8zp/6Vt27ZRv5tvvtniK6+8MvF15HEn8kryyxF1J3cdK005DKFRlm036PHzwqcQLV++3OL11lsv8ff2339/i994442o7dRTT7VY54NP09lwww0t9lVvlI6FVjQKIa4GQXWFb+j1S8c0hPhelXbfSkrXSEtHHDt2bNR24IEH1nmMcudXNVU2ypqmov06deoUtem9S48xcuTIqF/a0vvGTJcpgnKqQK2uTWk1Pf+8pGmMgwcPjtqK/ixRDk1ZSqsElHbd0XHLmnaWlf+7Y8aMsXjXXb/JgtCKgCHEqce/+MUvojafBls0aVsgpF3nNMVF0/SffvrpqN8OO+xgsY7P4sWLo35Tp061uG/fvlFb0nvBfyaphspfDUXPcfv27aM2nUdJqWohxNsF/OxnP4vasqar5eH+yUobAAAAAACAHOJLGwAAAAAAgBziSxsAAAAAAIAcapJ72qTR/RTGjRtnsd/L5OOPP7ZY91IIIVclv5oczRv0OZy6r4bmHt59991Rv48++shiLVXcvHnzqJ+OtS8z++qrr1qsY6972IQQ70/jyyKn5UWrM88802K/T4cvm5tXPi+3nHxNzQ8OIYRNN93UYn0vTJ48eY3/Flbl83JvueUWi7WseghxyeB7773XYv8+0PKKWeeDp7niet0dNGhQ1E+vHaibnytZ5045/Vq3bh216fhrqWKN66No8z5tfmT9v+r8+8EPfhC16VzUPZ/0Oafcv1W0saiUcudbGj3vSSWHQwhh7ty5dcZZjx1CscfVny+df3rPqc85yDonyjmv/nf0mfWKK66wuHPnzlG/fv36Wdy7d++oTffjrKY9wupDr5V+P8WTTz7ZYn3+958F9Lzr/n2e7oN0xx13RG2MT/l0Xg4ZMiRqS9o71e//d/zxx1usnzFDqPxzVENipQ0AAAAAAEAO8aUNAAAAAABADjX59Ci/HLR///4Wb7/99hb78muHHXaYxbpkP4R8LIFqqjTFwZcP3mWXXSzWcfNLEZNKRafxJYg17UmXJaYtH/bLF7Xk5rx58yzWkuQhxKV277nnnkyvtyh0SfLQoUOjNk1f0/M3f/78zMdnLmbnz5WWnz/00EOjtq5du1qcVg48K507Om9CiK8DZ599tsUvv/xy4jFQWUmpA55eH/v06ZPYT49RbnpU0WS9z6TRcfKlf7VN77O+PG0arqdrn15vjzrqKIt9Cr+m//hStnqN1feaTxnStOSijb2/jmVNidLz5eepXss01aUhyjXra5w1a5bFDz74YNRPP6vst99+Uds//vEPi1esWFHpl5hr+l7374WkUuH+3J5wwgkW77XXXhb7rRjS3jP6Ppk2bZrF9913X+rrR3Za5vu4446L2vR9oFsE+DSqSZMmWaxj1tSw0gYAAAAAACCH+NIGAAAAAAAgh/jSBgAAAAAAIIea/J42Wro2hBCuv/56i3VPjQkTJkT9pkyZYnHRcn3XJi2XPmDAgKitV69eFnfv3t3izz77LOq39957W7zvvvta3KZNm6ifjtsnn3wStem+M1paccaMGVG/hx56yGLNKw4hhKVLl1qs/y+/Z4eW7fRtTUW5c2DjjTe2uGfPnlGb5o2OHz/eYi3BiIaj7/sDDjggahs9erTFu+++u8V+f5KkPaV8TvDChQstHjlyZNR2ww031Pma2MOmstLy7rPOb71nbr755lGbjrmWzGyq17yGVO71VOffFltskXhMPf/Mo8aTtsde1jHX/ff0uuz3Fttqq60sHjx4cNSm11Td68jvvaLPMEV7zi33/5P2e3qNa8zzpePm912pra21uFWrVlGbfv7RfTuLNtYhrDr39P/u9/Pxe5gm/ftZZ51l8b333mtxjx49on46x/QZJoT42faCCy6weMmSJbuYtw0AAASaSURBVHW+BmSjpdUfe+wxi/Uzh6d7PN1+++1RW1E+d7DSBgAAAAAAIIf40gYAAAAAACCHmmR6lC4hHjFiRNTWsWNHi3WJoC/F3JRLfuWZnnNfinTixIl1xt6oUaMq+poqsaQ5zWuvvbbGxygCTTULIV7GOGzYMIuzlhz2P5MGUL45c+ZEP++///4Wd+vWzWJdLhxCCLvssovFWnb24YcfjvppWqqmboTAuK0N/rqW9Tqn5TPffvvtqK1Lly4W6303aSk66k+fbfySf13e/fnnn1usS/eRfzquy5cvt9jPUR3XPffcM2p78cUXLdZU8GXLllXsdVYDf86TSnunPUOmHa+c1+Hvn2PHjrVYP9+EEMK2225rsW4RUMTy3/7c6v/Rt+l4pY3JzJkzLT7kkEMs9unk+r544oknorYFCxbU2Q/14+9jRx11lMXbbbedxX4u6tYVmkZaxDkQAittAAAAAAAAcokvbQAAAAAAAHKopj7L+WpqanKxJfk222xj8bRp06K2Zs2aWfzVV19Z7HcD90u/865UKmVbn7kaeRnDKvVyqVTatRIHyss4Vjr1rCmkR1XTXMy6zLgJKtxcrAStHtW2bduoTat1vPfeexbr8uQQGvd9UrS5qEvEfTW+gw8+2OIxY8ZYPHv27KhfXq6T9dBk5mJahbZyzrtW0bztttuitvnz51v885//PGqbN2+exWkVjxpzG4CizUWl18UQ4rFPO8c6HuVeFzUtuV27dlHblltuafGbb75p8RpULmoyc7EhNfSWCg2tKc5FXz1v0qRJFu+xxx4W+/N/yy23WHz66acn9muC6pyLrLQBAAAAAADIIb60AQAAAAAAyCG+tAEAAAAAAMihJlnyu1evXhZvuOGGif20FKKWdgNQGZXOGy23VDEaBue/umjJ0rlz50ZtWoqW90XD0H1RXn311ajN/4zGV4k9StSMGTMs1n0bvCa4T1Gh+H1rksbe74XSvHlzi1euXJl6zCS6Z5gvi9yiRQuLv/jii0zHw+pxf2t8/hqn+9OqZcuWRT8PGTLE4moYN1baAAAAAAAA5BBf2gAAAAAAAORQk0mP0mWHtbW1FvslVboE8aCDDrK4MUsfAgBQJNWw9BhYW0iBqhz9vFCJ61bWY/h+K1asWOPXob83a9asqG3RokV1/i2gqZs+fbrFHTp0sHjYsGFRv08//bSxXlIusNIGAAAAAAAgh/jSBgAAAAAAIIf40gYAAAAAACCHauqTZ1lTU5O7pPZ11om35dGSpUVSKpVqVt9r9fI4hlXk5VKptGslDsQ4rj3MxUJgLhYAc7EQmIsFwFwsBOZiATAXC6HOuchKGwAAAAAAgBziSxsAAAAAAIAcqm/J70UhhDkN8ULKVdR0KKdjBY+VuzGsIoxj08cYFgPj2PQxhsXAODZ9jGExMI5NH2NYDHWOY732tAEAAAAAAEDjID0KAAAAAAAgh/jSBgAAAAAAIIf40gYAAAAAACCH+NIGAAAAAAAgh/jSBgAAAAAAIIf40gYAAAAAACCH+NIGAAAAAAAgh/jSBgAAAAAAIIf40gYAAAAAACCH/h9q/W1c6X8DjgAAAABJRU5ErkJggg==\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 9s\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": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0014, PSNR (dB) = 28.61\n",
      "MSE = 0.0036, PSNR (dB) = 24.40\n",
      "MSE = 0.0008, PSNR (dB) = 31.23\n",
      "MSE = 0.0022, PSNR (dB) = 26.61\n",
      "MSE = 0.0039, PSNR (dB) = 24.06\n",
      "MSE = 0.0006, PSNR (dB) = 32.18\n",
      "MSE = 0.0039, PSNR (dB) = 24.04\n",
      "MSE = 0.0055, PSNR (dB) = 22.57\n",
      "MSE = 0.0034, PSNR (dB) = 24.72\n",
      "MSE = 0.0033, PSNR (dB) = 24.83\n",
      "Wall time: 305 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": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9685\n",
      "0.9583\n",
      "0.9686\n",
      "0.9796\n",
      "0.9390\n",
      "0.9790\n",
      "0.9342\n",
      "0.9327\n",
      "0.9522\n",
      "0.9700\n",
      "Wall time: 3.46 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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9977\n",
      "0.9939\n",
      "0.9987\n",
      "0.9977\n",
      "0.9919\n",
      "0.9990\n",
      "0.9920\n",
      "0.9880\n",
      "0.9908\n",
      "0.9918\n",
      "Wall time: 11.2 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": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Simple Autoencoders in Keras - 784-64-784-tf1.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}