Deep Autoencoders in Keras - 784-256-128-64-32-16-8-4-8-16-32-64-128-256-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 = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_img = Input(shape=(784,))\n",
    "encoded = Dense(256, activation='relu')(input_img) # apply to previous layer = input_img\n",
    "encoded = Dense(128, activation='relu')(encoded)\n",
    "encoded = Dense(64, activation='relu')(encoded)\n",
    "encoded = Dense(32, activation='relu')(encoded)\n",
    "encoded = Dense(16, activation='relu')(encoded)\n",
    "encoded = Dense(8, activation='relu')(encoded)\n",
    "encoded = Dense(encoding_dim, activation='relu')(encoded) # the bottleneck = encoding_dim\n",
    "\n",
    "decoded = Dense(8, activation='relu')(encoded)\n",
    "decoded = Dense(16, activation='relu')(decoded)\n",
    "decoded = Dense(32, activation='relu')(decoded)\n",
    "decoded = Dense(64, activation='relu')(decoded)\n",
    "decoded = Dense(128, activation='relu')(decoded)\n",
    "decoded = Dense(256, activation='relu')(decoded)\n",
    "decoded = Dense(784, activation='sigmoid')(decoded)\n",
    "\n",
    "autoencoder = Model(input_img, decoded)\n",
    "autoencoder.run_eagerly = True # plaidml / keras not running eager mode by default\n",
    "\n",
    "encoder = Model(input_img, encoded)\n",
    "\n",
    "encoded_input = Input(shape=(encoding_dim,))\n",
    "\n",
    "decoder_layer = autoencoder.layers[-7](encoded_input)\n",
    "decoder_layer = autoencoder.layers[-6](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-5](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-4](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-3](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-2](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-1](decoder_layer) # extract last layer\n",
    "\n",
    "decoder = Model(encoded_input, decoder_layer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1595257946\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 19:22 - loss: 0.6930WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (4.808486). Check your callbacks.\n",
      "235/235 [==============================] - 92s 390ms/step - loss: 0.2806 - val_loss: 0.2489\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 64s 271ms/step - loss: 0.2379 - val_loss: 0.2199\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 73s 312ms/step - loss: 0.2129 - val_loss: 0.2062\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 80s 339ms/step - loss: 0.2022 - val_loss: 0.1972\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 76s 323ms/step - loss: 0.1944 - val_loss: 0.1914\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 83s 354ms/step - loss: 0.1902 - val_loss: 0.1885\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 84s 358ms/step - loss: 0.1861 - val_loss: 0.1839\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 88s 374ms/step - loss: 0.1830 - val_loss: 0.1824\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 84s 358ms/step - loss: 0.1812 - val_loss: 0.1817\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 84s 358ms/step - loss: 0.1805 - val_loss: 0.1824\n",
      "Epoch 11/50\n",
      "235/235 [==============================] - 89s 378ms/step - loss: 0.1790 - val_loss: 0.1785\n",
      "Epoch 12/50\n",
      "235/235 [==============================] - 83s 355ms/step - loss: 0.1773 - val_loss: 0.1772\n",
      "Epoch 13/50\n",
      "235/235 [==============================] - 91s 387ms/step - loss: 0.1766 - val_loss: 0.1769\n",
      "Epoch 14/50\n",
      "235/235 [==============================] - 86s 364ms/step - loss: 0.1760 - val_loss: 0.1756\n",
      "Epoch 15/50\n",
      "235/235 [==============================] - 72s 307ms/step - loss: 0.1756 - val_loss: 0.1759\n",
      "Epoch 16/50\n",
      "235/235 [==============================] - 70s 300ms/step - loss: 0.1751 - val_loss: 0.1752\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 69s 292ms/step - loss: 0.1744 - val_loss: 0.1756\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 66s 282ms/step - loss: 0.1739 - val_loss: 0.1737\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 64s 271ms/step - loss: 0.1731 - val_loss: 0.1759\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 57s 241ms/step - loss: 0.1731 - val_loss: 0.1739\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 60s 254ms/step - loss: 0.1721 - val_loss: 0.1720\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 70s 299ms/step - loss: 0.1713 - val_loss: 0.1727\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 57s 244ms/step - loss: 0.1713 - val_loss: 0.1714\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 53s 226ms/step - loss: 0.1710 - val_loss: 0.1729\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 60s 255ms/step - loss: 0.1712 - val_loss: 0.1729\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 51s 215ms/step - loss: 0.1716 - val_loss: 0.1743\n",
      "Epoch 35/50\n",
      "235/235 [==============================] - 59s 252ms/step - loss: 0.1748 - val_loss: 0.1776\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 52s 220ms/step - loss: 0.1762 - val_loss: 0.1750\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 57s 242ms/step - loss: 0.1707 - val_loss: 0.1702\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 57s 241ms/step - loss: 0.1688 - val_loss: 0.1692\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 49s 207ms/step - loss: 0.1683 - val_loss: 0.1682\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 39s 167ms/step - loss: 0.1628 - val_loss: 0.1637\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 40s 171ms/step - loss: 0.1616 - val_loss: 0.1622\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 39s 167ms/step - loss: 0.1606 - val_loss: 0.1614\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 48s 206ms/step - loss: 0.1600 - val_loss: 0.1611\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 34s 145ms/step - loss: 0.1593 - val_loss: 0.1602\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 41s 173ms/step - loss: 0.1591 - val_loss: 0.1607\n",
      "1595257946\n",
      "Wall time: 53min 9s\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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoded_imgs = encoder.predict(x_test)\n",
    "decoded_imgs = decoder.predict(encoded_imgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": "iVBORw0KGgoAAAANSUhEUgAABG0AAAGMCAYAAABgTl9HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAANgUlEQVR4nO3bsWqX9x7H8SeH0ppKSIMpBEI4WWwvoJkKXoGDi7uTg1269R7qBehgl44uQtdC1045Q4cOpRlaoSWixBpBpQ5PL+Ao/MP5ef5v/75e88OHBx6+Bt8ka/M8TwAAAAC0/GvZLwAAAADAfxNtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAILeO8vD586dmzc2Nt7Uu/AaT58+nV68eLE2YuvChQvz3t7eiCnO6Keffno0z/PHI7Y++uijeWdnZ8QUZ3B8fDz99ddfQ25xa2tr3t3dHTHFGf3888/DbvH999+f19fXR0xxBs+fP5/+/vvvIbfo39Pl+eWXX4bd4ocffjhvbm6OmOIMnjx5Mj179mzILZ4/f37e2toaMcUZ/fHHH8NucXt7e97f3x8xxRn89ttv06NHj4bc4vr6uv/zL8nDhw9feYtnijYbGxvTlStXxr0VC/nuu++Gbe3t7U0//PDDsD0Wt729/fuorZ2dnenOnTuj5ljQ9evXh23t7u5O9+7dG7bH4j799NNht7i+vj59/vnno+ZY0I8//jhsa2dnZ/rmm2+G7bG4S5cuDbvFzc3N6dq1a6PmWNC33347bGtra2v68ssvh+2xuK+++mrYLe7v70+Hh4ej5ljQwcHBsK2NjY3p6tWrw/ZY3O3bt195i/48CgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgKD3zvLw48ePp3v37r2pd+E1Tk9Ph23dv39/unHjxrA9luPo6Gi6cuXKsl/jnTPyFs+dOzd98sknw/ZYjpcvX05//vnnsl/jnfPy5cthW7/++ut0+fLlYXssx8nJyXT37t1lv8Y75+TkZNjWgwcPpps3bw7bYzmOj4+nr7/+etmv8c45Pj4etnV6ejp9//33w/b43/lNGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAICgtXmeF394bW3xhxlqnue1ETu+4VL9Z57ngxFDvuPyuMWV4BZXgFtcCW5xBbjFleAWV4BbXAmvvEW/aQMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABL13loc/++yz6fDw8E29C69xcHAwbMs3XJ61tbVhW77jcoy8xQ8++GDa29sbtsfijo6Ohm25xeXwc3E1+Ln49nOLq2HkLV68eHG6devWsD0W88UXXwzbcovL87pb9Js2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAEDQ2jzPiz+8tvZwmqbf39zr8Br/nuf54xFDvuFS+Y5vP99wNfiObz/fcDX4jm8/33A1+I5vP99wNbzyO54p2gAAAADw/+HPowAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACC/gFpgTMSLVjEkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2d29dWVfXHJx1MNH+JAgKChEcQFEVOQgmYqKWji0wd4WjUTY7G6C/oT+i2ixyjBmkXWRYdDMvykBhHBeLkATnKSeAVRcDKjvwuGqy+8wt7+/LyHvaz38/naj6sxX6ed6295lp7j/mdc8jJkycDAAAAAAAAAACaxUcG+gcAAAAAAAAAAMDp8NIGAAAAAAAAAKCB8NIGAAAAAAAAAKCB8NIGAAAAAAAAAKCB8NIGAAAAAAAAAKCB8NIGAAAAAAAAAKCBfOxsOg8ZMoT64APEyZMnh/TGdZjDAeXIyZMnR/TGhZjHgYO12ApYiy2AtdgKWIstgLXYCliLLYC12ArOuBaJtAHoP/YM9A8AgIhgLQI0BdYiQDNgLQI0gzOuxbOKtAHoC4YMyS+FT54888td76ef//Of//T+DwMAAAAAAAAYQIi0AQAAAAAAAABoILy0AQAAAAAAAABoILy0AQAAAAAAAABoIOS0gV7lIx+pfg/4iU98otjnnXdesS+88MLUT/PTvP/++8X+xz/+kfp99KMfLfY///nP1Pavf/2rm78YAAAAAAAAoJkQaQMAAAAAAAAA0EB4aQMAAAAAAAAA0ECQR8E58fGPfzx9VtnTJZdcktqmTp1a7OnTpxd73Lhxqd///d//FVvlVu+8807q9/bbbxd72bJlqW3jxo3FPnbsWLFdNqXlxatKjQNAXosug1Tpo6PSRV1/KoOE5jFkyJD0Gf8IAACdiu9p+rnufFPX9rGP/e8xWvdITd/gn/13VJ2RPCWEp4GAwQeRNgAAAAAAAAAADYSXNgAAAAAAAAAADYSXNgAAAAAAAAAADYScNnDWqDbz/PPPT22aj+byyy9PbRMmTDijPWvWrNTvyiuvLLZqP0+cOJH6vfnmm8X2/Bia/0Z1oX/5y19Sv3//+98B/YvreT/5yU8W+6677kptWvJ95cqVxfZ5JD9K7+A67EsvvbTYkydPLva0adNSvzlz5hT74osvTm2bNm0q9tKlS4u9ZcuW1O/o0aPFZl12j+5q8hUdW89To/+vTpOvucz8GupvfV3qd7Nmewf3p8OHDy/2ww8/nNouvPDCYv/whz8s9t69e1M/zZ1ALqOBoWpte1uVHVG/1snndzo6foxJO6jb0zQHp7Z5jj7tN3To0NSm5x39f3qu9ba//vWvqe3IkSPFPnToULH1/BvB/tnfuD/Ve6ku519fnl+JtAEAAAAAAAAAaCC8tAEAAAAAAAAAaCDIo+CsqQvH1fAxL6+t8iYNATx48GDqd8EFFxR7xIgRxf7Upz6V+qlcY/fu3antD3/4Q7EJc202KsG57bbbUpvO3csvv9xvv2kwoWv2sssuS23z5s0r9gMPPFBslUNF5BBhX29jx44ttvqL48ePp34aCow86n/omHl4t8pd1Ha/rOg4e0lRLV+qdkQOC1cZrP8mnbuurq7UprJGleAQ6t1zfK51z/zCF76Q2lQ69eyzzxbb5VHsmWembl1Vhc7XSRXrxlkliGpHZFm6rktfR3/729+K7ZKMv//972f8f526FnXM1Xf52NXJZVzufwodq4g8bzpedSWlfU/Tz3pW9nMza/HsqZIK18mjhg0bVmzd3yKyT/W2iy66qNgunVL03tL0DRF5L9Q1+8EHH1Rer23UrR1dA3WypLpn06r/U3cN993qS1xCp9fQNex+91zXM5E2AAAAAAAAAAANhJc2AAAAAAAAAAANpF/kUd3NdF8XRqp0N7xIwxbr/k93r1cXbtXTa3YiOq4eXq+hfZoRPSJi48aNxdax1BDFiIhjx44VW+UaKoeKyGGJXqnqkksuKbbeV3VhcTAwqLRm0qRJqW3nzp3F1lDRTg3hbgq6/lR26PK0b3zjG8WeOnVqsb1ClK4xDfWNiBg5cmSxVa6hviIi+wvkGv+jzldq2Lb6Su+nodkqU/XqFIrLozS0XP2ryw/qrq8yA5cBQM/wtTF69OhiX3vttalN7yXdu33NDmbq5IjdrTCjbS6ZUBmjVphRWXhE9st61nH0N3lVRZWNu4S80+WJfh5XP6Tj6mOnEsGrr746telepfPp/lS/S2VOLqNS6b+eZSIi9u/fX+z33nuv2L6e6yqADWa6K1XUfczX2JgxY4p9/fXXF/uqq65K/VTi7deo2sd8/1R/q5VvI/Izj0qn6qR9nbJm6+RGOkbuT6vavJ/K1VTi5vOkaRj8Gno+0nHViqYREW+//Xax/flW29wP9yZE2gAAAAAAAAAANBBe2gAAAAAAAAAANBBe2gAAAAAAAAAANJA+yWnjujXVi6lO10s4a54E1aRG5HwWqut0zbFqAFVX5qX89PquPdTvUi2rl+xTDduuXbtSm2rhuptbp1Ooy2mj+Qy8XJ2WPlPtoWtCN2zYUGzVIX7xi19M/WbNmlVs136qNrmqFCc0gxkzZhR74sSJqW3VqlXFHkzlD/sa9Xmq337ooYdSvylTphRbfaFrfQ8cOFBs1WdHZM2x5nO4++67U7/t27cX+/HHH09tnitgMKH+y/fMT3/608W+5ppriu15HHTt7Nmzp9i+H+l+qnMVkfPY6P7pe4D+XteOt20vbAK+p916663F9vK0mi+h7j4YbFTlsfFcJpqfxs+oesYcNWpUsa+44orUT/NjqB/2M6quHc8hpr9Rc4N5rgw9A6uPjuj8s5D/fh0v9VWet2bu3LnFnj17dmpTX6a+qq5cuubt8tyKes7V9RYR8cwzzxR7+fLlxT58+HDqh8/8cOpypeg69Vw106ZNK7Y+T4wbNy71073QnwP1nunq6iq2n1f1fvX92fOvnMJzjXViHhv1Vf6cVudP1ffq85zmbIuIuPLKK4utvtXHWPdCP5coei/52eatt94q9tatW1Pb6tWri61rva5EeU8g0gYAAAAAAAAAoIHw0gYAAAAAAAAAoIH0mjxKQ4C8xKGW2LvpppuK7WGLGrrmoU0anqbhS1rmNCKHHmkIY115RpfneCnaU3j50jfeeKPYv/71r1OblrfW0Mq2hTd6uF5dGUmdA5WWeQiaSi80vFTD/yMiJkyYUGyXT+g1ta1TwgvbjMsR58+fX2wPR//5z39ebOau53iIppaKvv/++4ut/jkiz4eWOFy7dm3qt27dumKrRDIih6/eeeedxXbf/bWvfa3YL7zwQmobzFIOnYPx48entjlz5hR70qRJxfZylLpXafi9hnNH5L3VJRkaoqxzV1fatE6CMdjmsa/wkHNdYy4f1xBun/vBTFVZWh8/DeHXc21ElkFNnjy52C7JUCmEzoGG3kfkNezyD5UI6G/ys6zK4fwM7NfsNPz36zpQKYTL09RPuixFfdLevXuLrdLdiDxXeuZVfxyRfaZLv9UPv/7668V+9913Uz+9vvvMtvvQ7pb19nWqz48qdVq4cGHqp1JSlTT6OVQlUS6V0/nS0u1+XtXf6L9X51HXsF+jU+Zb501tH1ddfy7l1XlTH6ryxogsiVIfoL4vIktFfe/Ts46m5HCZnMrT/T2BnrE8XUpv0tleGwAAAAAAAACgpfDSBgAAAAAAAACggfDSBgAAAAAAAACggfRJThsvS6pl8DQviWvTVFOoJfsisvatrly351M4heu+VaPo/0c/629yVDepOR0iIl577bViuwayzajm0vWXqs1V3aZrD1UrqHPtOkTNi+Ol8XQO9d7sFE1om3GN+cyZM4vteaN27tzZL7+p7bj/mz59erE///nPF9vzmKhGW/PYPProo6mflpN1Hba2aR4qL92ouQZuu+221KYlwD1vQ9twHb/ufarBj4j47Gc/W2wtGXzo0KHU75VXXim2lgVWnXdEzhPg46wlTHVf9HtL/bT7Wy+XCueOln2PiLj++uuL7WtR1y1z8eHU5U3RXAoREVOnTi32lClTiu25ZLRUrOY+3L17d+V3e54TPXvq73D/XZe7Uee/DecifRbQvFpeSl37+blR5+PVV18ttuacicg5bfSMqjkYI3JODM3tFpGfi9Q+ePBg6leXJ7IN81aXt6auTfPC+Bxrzr6bb7652Lfcckvqd9111xW77plEzzD63BERsWPHjjNew3+Trtm6vDj6/OP3Zyei8+T3r86vP2trDqgbb7yx2J4jTM8lmktmw4YNqZ/Ok68bzRV4xx13FNvzLmq+G1/r+rfo39nba5RIGwAAAAAAAACABsJLGwAAAAAAAACABtJr8qiqkmURWfKgYWdaxjWiPjRbS35reJqHd2tpPg0r91A1DUdzeZRKszTMzkuSaVhkXQm3NoQw9gY6vxq+6iW/NWxX511D0yLy/Hq4oc6Vhid7mTZCxPsHDR3UkrQReb39+Mc/Tm1eyh26j465S5EWLVpUbA83VTZv3lzsRx55pNibNm1K/XSeXFKg4asvv/xysVWiFZHDTT2MecmSJcVuuzzK9xKVE6ocKiLPncorPJxfJVG6H3v4tc6d+0aVhmjYsIcJq3TA16+HR0PP0Hl68MEHU5vufR7Kr+uZve9/dPe8ptJ/91/qs4YOHVpsL/+qPlAlOMePH0/99HzppWf1PKvz7edh9ZXqhyPady6tkpL4uOoziPsjvYZKYvbt25f6VaU8cFmNSkVdHqXzq/uznnkjIo4dO3bG72oLfh9WlaL3f9d9UqVHEXlf1DLsKs+OyPOlzwbbtm1L/fS50udD7y+9f/z3qnTR90W9ht5bbfDROr8ud9Ozg68PlUfpGchTKKj0WyVRLo/S+XUZqZ5n9L5SPx6R/YX7Dn2O7UtZG5E2AAAAAAAAAAANhJc2AAAAAAAAAAANpNfkUYqHj2l4lFYg8aoTVZWfInK4kVY00RDGiBwOrCGkLsE5evRosb3a1e23317s+fPnF9srVWm4lUu9NBS1bWGoddSFGWs4mYYAehihhvtqZQwPn9PQtbox1vvAv6sN4YedgFax8Ao4Oo9PPPFEakNO0XNUFjht2rTU9rnPfe6M/bxa12OPPVbs9evXF9vDw3WePARWfa1W1fN9QqWoWi0gIq/htlfjcynvggULiq1VabyvyqNUBhOR90zdC319qX/0qjcaQqzVFryfzo+HgQ+mvbAvUWnFvHnzUpueq1R+E3G6fANOR+9RlWBHZD+q1aIicgWg/fv3F3v16tWp35YtW4p9+PDhyt+hIfxeIUzlH9rPq0zp+vMzcKevRT+7qd/R6nkqbYnI1dVGjhyZ2lSuUVd1VJ8F1Ad7Px1zlyrq/1NJlO+fddWj2oiOoe5HLhtWSZTLB/WMqXuVy1aOHDlyRtvvGZUv+XNqVXUkl9bovLpkT6/Zl5WHBhr3p+oz/flOJYO6jnxPW7lyZbHVt+qZJyLvi8OHD09teq5SX+tSRb2mnmsj8vpWWSrVowAAAAAAAAAABgG8tAEAAAAAAAAAaCC8tAEAAAAAAAAAaCB9UvLbtbOqA1MtqOvuVW/t5Qm17eDBg8X2EodVpbbqysSqri4iayVVL+zlojWHwI4dO1Kbj8FgpE7fq/eBl4xVTeGsWbOKffXVV1d+V11ZRNUyul4Y+gfVjs+cOTO16dp0LTH0HM3V9aUvfSm1aZl11VMvXbo09VuxYkWx1de6tr6urKOWaFTdr2uCFdc+e/6ztqHafddbqz7f2zR3w5o1a4qtJb4jqnOs+Vzp7/DcOmPHji223j++L2pugLp9t84Xt03L39uo/v/aa69Nbbo2n3zyydTm5yr4L1VrwvMd3nDDDcX28sHqo3T9aR6viOqziubyi8h5Ojx/jn63+lfPaaM5c3y9VeXicJq6Fv13qa/RvUrzC/nnm266KbVVlVJXfxeRx1LPNpdeemnl7/VnE895cqbfHjE48tgoOrZq+3jpuE+aNCm16TOEniU896iefTQnkpeE1n3Ry4tX7XFV8xtxeq43zQtXl0upU6jKbaq52CLyPGlZ74jsT/fu3VtszWETEfHGG28UW3MFeQ4kXc8333xzaps+fXqxx4wZU2w/26h/1ed//269X9zvnuucEmkDAAAAAAAAANBAeGkDAAAAAAAAANBA+qTkd13YooajeQk8bfNyXdqmMhsv+9fd0CMNg/SQc5dvnMLDFp999tlie6j/YAtp7A4qidIQQw27j8iSKA391rLRETks2Etn6j2ic+0hcz0JRSSs/+y54447iu1h5VoS1UshQvfxcvZaonnGjBmpTe9hXTvPP/986qd+rbulKL1N/5/6a1+L2s/lUG2XNerf6zJQlUn42GposJZr91D8qnBd9ckRucSl3j8Red2qXGDfvn2pn+7dPo/qw/Ve8P1S2/Cp/0XXy3333Vdsnyc9p6i8MeL08xL8F10Tes+6dF7Pih7qr/IKtb2fzpeuBy3jHZHPoX4m1VD/rq6uYvu5WcP7Xdah94L6B/cddRLHJqE+RH+zj4nKGrxNZSpaGtzHTvupTMclHqNGjSq2pwHQ+0xt98ltp06iqxIjnwOVMLkPVFljlTTYP+s69VLPI0aMKLZLm3S+9Pe6vFhlee6HtYy1rr82nHt0fPxZW+WEfn5Vv6ZnGy+53l2p4nXXXVfs+++/P7VpyW/1I+rHI/I8uey1KiUK8igAAAAAAAAAgEEAL20AAAAAAAAAABpIv8ThaTiQhpZp+FNEDiPyEE393BPpkYdeaWidZ5LWjNb6vWvXrk39XnjhhcrfO1hDuutC7zVUVLO9T548OfXTSgk6T14ZQcPTtm3blto0/FCzx3uoclVI7dlQlSl9sKPh/Bre7TK3X/7yl8X20FPoPh6GqZIoD/PXe13XjlfvUr/W03u7ShbpIcjazyusdEqI/tmg86U+SsN4I7IUwsehqhqQV7/Qcdf152tR7xPfF9Vna+i3hxBrm4eI6/2k4cS+7vGjp6PrZc6cOcV2CZqGknuFC/hwdC26JEPH2uWderbQ6l4ubdJrqGTG16zKJF1WoNfQNeXVTvSaft5Wv6JruFMlylXnMPctuj48RYOe/VUyMX78+NRP50Olol4BTO8f993q//V76yoE+jXa4Cf92UzXle4fdRUlfY51X9Q0ClrZMCLL43TtuE/VdeTnFk3vMHr06GL7vqjnG62G5L9XaYM8SufTnwlVyqsVoh0dS5Wq+We9vq6piIi5c+cWe9q0aalN51T9pJ+H//znPxfbZeF6nunLdUmkDQAAAAAAAABAA+GlDQAAAAAAAABAA+GlDQAAAAAAAABAA+n32nKaa8T1t4rnrelu3hDVAKqtJfoist7tnnvuSW2qUd21a1exFy9enPqpHrYN2tLeQMdZc9hERNxwww3F1jw2rj1U/bBqVV0HquVuveS66ry1/J+XZdOcDq4r1XtQr1eXv6iq7NtgRPXI06dPL7bryJcsWVJs1lHP8RwLmktBS1FGZN+7cuXKYnsJ1N7IH6ZadC1rq37Wv0t1/BHtXFe6P3nuF0V12t5Pc2fcfvvtxX7zzTdTP52DurLe6r8nTpyY2saMGVNsvdcOHjyY+qn+vK7Eqs53T+6zwYaWM7322muL7fvRj370o2JrPgeopip3hOfW0nwvnrNCfayuS8/BoP3UV/pa0XXq/rvqzOE5O/T85H5Zr1F1bvbPnbI/qz/x8sp6/lizZk3lNXR/8vLBmrtE+3m+Ez1TdjcHkuZ0jMh5NdqS868u96Xe6zru/jyhY+33pea71Pn35wT9rGeiulxv48aNq/wdus96Xhz9rrpnCM1b5Pdup6w/RefafaauRX1O8/9Xlf8pIp9ZdC48N6Dm5KvLbXr48OFiP/nkk6nf6tWri+33kq7NuvPMufpTIm0AAAAAAAAAABoIL20AAAAAAAAAABpIv8ujNGyoLnS6p1SFfGppxYiIBQsWFNtLm+pvXLZsWbHXrVuX+nmI22CkrjSehqNF5JBuDV3zkFItuaYhaCpVi8jhjF4OUK+p4ateVlPD8zykTcOOT5w4UWwviall6/wabSxV3F00jFTLlz7zzDOpn5cihZ7ha1HDSD3kXe9hLYHaG+HX/ju0hKnK5DwcVkOBtQx5RDt9bZU04sCBA6mfhnr7mGm50XvvvbfYLnPTsdV7wcPA1c/5nqkhxSq78bLSKmP1cGj1qb1RTr7NuJxl/vz5xVbJjUtiVG6K7Ozs0XvU1+KKFSuK7ZIZnROVePj5Ru/1OrnVFVdcUWyXkOg63bx5c7FfeeWV1G/Pnj3Fdqmc+no9S7VBklF33+uYL1++PLWpNELl/CpNjMhjqedLl/Hqd11zzTWpTcuD61xff/31lf30rBnRuedLlYr5HqRn9Lqy97rGfL6rzjS+L+rzhV7Df5Pie7BLF0/h60j3O12/EXleu5sKpMlUSYB0fTn+HKD3gd4v/vylz366xlwCpZ99XHUPffrpp4v94osvpn66H/j86mdvU851Tom0AQAAAAAAAABoILy0AQAAAAAAAABoILy0AQAAAAAAAABoIP2e00bpa72e6hK9fOndd99dbNcovvXWW8V+4okniu16yDZoD3uC6hUvuuii1DZlypRiawlab1Mdout0tcym5rHxvBb63Z6zQ0vjaptfQ7/LdaaqnVQto+fs0JJzrjHuVM1xT/AcDA888ECxVZOvmtGIwTVGfYmPv2rh3T+pL1M9r89Fd/2azq/7U81jM2/evGJ7rgctHb1q1arU1saS3zq26gM1R0VExK9+9atie94Lz5NwCs+BoXr9Kv8aEfHOO+8U+6677kptqtPWXGCe623fvn3Fdv05Zb67j2vyFy1aVGzV+GvOo4icxwS6h96LWqZZ10NExEsvvVTsQ4cOpTYt26w5Sjy/gfqyd999t9ieb0+v4WtFzyPPPfdcsT2njfoL96Hq6+vKDHc6dbkKPWeRngHVN/qepvlV9KzpY6znzTvvvDO16flV92rPaaPlptXvRuS56qRnED2T+/zoeI4ZM6bYXmpbz92ei09z1eicau6miHzfa56iuvLT+r0R+TlT+/m9oPea7pEReZ9sQ643/d11OYV0zPX8F1FdPt3Hf/z48cXWvKmaayoiz5Pnz3nqqaeK/dhjjxXb91Kdm7rzS1Vu3d6ASBsAAAAAAAAAgAbCSxsAAAAAAAAAgAYyoPKovkAlApdddlmxv/KVr6R+V111VbE9tO5nP/tZsV999dVi91Q60DZ0jD2kV0OEteRwRMTo0aOLrWPuIcgaKqihkhpCGpHD3Vxqoeg1PKRNQ/C8JKbOfV25PpUbeGhxVem7NuKlD1WCqOGIXmYTegeXR1WVkYzIvkzvbQ/l9GuewuU36ge0VGpElnVoyL/702XLlhV7w4YNqa1tIfsR1SHEHrr7+9//vtguf9CS31py2NGw8DfffLPYGkYekcsYL1y4MLXpHOzYsaPYW7ZsSf3UH7Jnnh26/lwOoLIJlfBoie+I0yXAcHboPet7vY67l69VWeOwYcO6dX09wyxYsCD10zLTLrVYu3ZtsVWydTbymSp5f9vXqP597p/Ud6m0zM+oOs56Dd8/dZ/0EvEq69D0DS7N0X5vvPFGavNnl05B70sfMz3njxo1qti+v+n/831MJTR6LtX1G5Flpip501LjEflMc+utt6Y2/V36DLF///7UT8+9WpLc/1+nSt6UqjXmz196/7p0TedUn+F8Hek9ovOkZ6OI7EP/+Mc/prbvfve7xVZZpPuHOklUlQyqt+eQSBsAAAAAAAAAgAbCSxsAAAAAAAAAgAbSOnmUhsLdf//9xb7nnntSPw3B85DDX/ziF8XW8NhODVXrbTQMzCUxKlPy6hfapmGJXo1J/5/KLjRELiKH1rk8Sq+v4eKeeVzlHx6ep99XF1Krn+tCkNuOhy3q3G3fvr3YXnUDegeXRezevbvYs2bNSm2aWV9lF3v37k39NGxX14pLFfUaDz/8cGrTcGJdl9u2bUv9NGu/VlVpK1U+xaVgGtLt46JyKQ3Fd1mbXlO/S+cjIvtiX8/6e3UNe0WINlS/GCh0b505c2Zq0z1Ox1xlhXDuaAh8XTUgPy/o3Kl0yteY9lMpjMpgvM2rxv3pT38qtlb/q6uyx1o8HR+TKmmKyyJ0Duuqw6iv9UpV69evL7bKhl1ap1IdnXf/7k6a3zo5ns6Bjruf3fVv11QYEfnsrmvRz57aT9M5aBWiiIhJkyYV22WrOsdaeVgljBFZ8u17ZhskUVXUVYzUOfRzj55h1IdqapOIiNmzZxf7xhtvLLZL+Ddu3FjsxYsXpzZ9B9Ab5xeqRwEAAAAAAAAADDJ4aQMAAAAAAAAA0EB4aQMAAAAAAAAA0EA6PqeNa/e1dN43v/nNYrvmUUv2aQ6biJx/o42lZs8V1eh5yUHVamrJ4YjqXDU+h6oprCo5HJHzDWmpRr+G6ry9rLf+Py+XqaXf9uzZc8bvjch/c11JuDai94JqSyOyplTLArdNs9sU6srCfvnLX05tuv6++tWvFttzQ6kGXHMsTJ48OfX7zGc+U2zPzaC6cdV8P/LII6nfunXriu15o9pIX5bZ7a6O2vOJ6Vx53qIq3b3n9mB99xzN2/D1r389tWl5aC1BTInv/qMna7YuV4OuvxkzZqR+msdB89ZERGzdurXY6vdZe+dGVR4bP9epz9NzjvtdXZueT0WvoWvbS37fcsstxX7yySdTm55ZO+nsWZcjUs/Tmo/Gn+H0s+df03OMPqP4M8kFF1xQ7DFjxhRb89tE5HXq+Wg0d+DKlSuL7bnGurq6zvibIvLcDaY1XPe36lrSHIxz585N/RYuXFhsPbMcPHgw9VuyZEmxX3755dTW2z60L3NNEWkDAAAAAAAAANBAeGkDAAAAAAAAANBAOl4e5eH83/rWt4qtIW4eovTiiy8W+6c//Wlq0xC6wRSq1l10TLwErUrLVO4QETF8+PBia9k2Lxuu4ZIatuZlL3fs2FFslTJFRBw/frzYOp8uj9I2D1nct29fsTUc3eUAgzk8WcMAvUyizo+WqqwrSwo9x8PwNTzX18eUKVOKPXXq1GKrbCoiryMNUQx3fr0AAAdmSURBVPVwZF3DLtdQaeH3v//9Yv/kJz9J/XxtQs/pqR+66KKLiu2h/uqXNZwbCfG5oeOs4foTJkxI/dRvrlq1qti+B0OzqFuLI0eOLPbll1+e2nSN7d27N7WpLGawnTn6Cx3XOh+n/VzOX1c+XqWomjpAz8kR2Q+MGDEitWmp4k4t/+1jq/IjLZPt+5H+jf4cqOtKzzReNvyDDz4otu59+u8RWQK1Zs2a1KbnXG3buXNn6qfPDS5l66T56it8fnVO58yZU+wHH3ww9dO51mfEp59+OvVTeZTL5Ppy/Ovu255ApA0AAAAAAAAAQAPhpQ0AAAAAAAAAQAPhpQ0AAAAAAAAAQAPpyJw2qhvVUrMREffee2+xtWSil9v7wQ9+UGzNXRKBRv/D0NwGnodi27ZtlW1a7lfLW7rOVDX6WvJv//79qZ/Om5frrioZXFdqz3Nx6DXU9n6D+X7RUoiu3dRcJq+++mqxO6k0ZSfhWlnNg/Cd73wntelnLdGtuW78mmr7OtL1t3HjxtS2ePHiYj///PPFdl0x9D++ZrVkpvs1zcmgftl9Lfr8njN27Nhi+zhqXjXNV+U51qDZ6LlUc/tpXpOIvN48P8axY8eKzXrre3yM1TfW5bTRuVbfGhExatSoYmuZbz1TReSS4pp3xa/v59JOwc+Dei7QM77/fQcOHCi25jWJyH5Uc0X5HOg60jPMkSNHUj99dvFS0lVl3TUfoPfjDHw6ep9HRIwePbrYixYtKva4ceNSP/WTK1asKPb3vve91E+fK/t6/PvSJxNpAwAAAAAAAADQQHhpAwAAAAAAAADQQDpSHnXeeecV+7777kttXor2FJs3b06fVa7RqWGFTcDLGGqI4dGjR1Pb1q1bi63l1xydDw1D7Y2QNpcD1EHY8ZmpKlH72muvpX4aUqpymrqymNB76Jj/7ne/q2x76KGHij179uzUT8suarivyiAjIp566qliL126NLWprLFKtggDg69FXdsqb4zIofm6tt2ndmrp2YFC52Do0KHFdn+qofwqfWSMm42vDy31rHJUPddG5PK1u3btSm16RmK99T0+rjrmaqtcKSJL3lzWoVIdvYbOe0Quge3X13vJy1R3Cn7+03tbzxz+9+ke5CkW1Kf6HqdUpUDoaT89V7m8WP9O1ul/0fvZ5aHz588v9qRJk4rt/lTLsT/++OPF9vNLW57zibQBAAAAAAAAAGggvLQBAAAAAAAAAGggHSOP0pCoK6+8sthePUpD4bTawnPPPZf6adgdoWq9h46lh6M1ITyNuT53dAy1QtiaNWtSP80GrxUBkEP1D1XzFBHx29/+ttjLly8v9pgxY1I/lb+pFFJDtiNytQWvCsWaaxZ1ElGtjLF27drUpveCVrNhPZ8bGka/ZcuWYn/7299O/XSctXpXE/ZVqMbXm0o51Fd2dXWlfio1dzlqlTzKvwvf2zt0V1bvlZ9UAuWpBF555ZViX3zxxcV2GZVK4+r21rbI5NQfqiTKK1a6VKw7+Lj0xjh1dw46eU56C19H+rw+bNiw1DZx4sRi67OEVxHWs6xWj/Kqin09/v01v0TaAAAAAAAAAAA0EF7aAAAAAAAAAAA0EF7aAAAAAAAAAAA0kI7JaaNa0RkzZhTbdXBaLm/Tpk3FXrduXerXqeXxAJqEausPHDgwgL8EzgYtW/nuu++e0YZ2otprz7OwefPmYmtOuIhcknPfvn3F7m/teJvRnBUbN24cwF8CfYX63u3btxd7/fr1qZ+eWTXXUUTO78F663vqcnFobhUvPa3nI80XFpHzwOlcjxw5MvXT3Eavv/56atPnnTbeB/o3+d/XxFxqbckrNBB4KXW979VP7tixI/VbsmRJsXVNecn1tkCkDQAAAAAAAABAA+GlDQAAAAAAAABAA+kYeZRy/vnnF/ull15KbUOHDi32b37zm2Jv3bo19dPQKcLYAABgsOEhxCrP2b17d2pTSYDKM5oYpg7QFHx9nDhxothaotbLeh87dqzYLlVU2Q3rr+/xMdbPKuvojbQL6mcj6iVC0CyYn3p8fHTtdHV1pbZHH3202EuXLi32oUOHUr/33nuv2OoX2zoXRNoAAAAAAAAAADQQXtoAAAAAAAAAADQQXtoAAAAAAAAAADSQIWej+xoyZMiAicS0rJ7mrdFS4BG5hKnq7l2732l6t5MnTw758F4fzkDOIcT6kydPTu+NCzGPAwdrsRWwFlsAa7EVsBZbAGuxFbAWWwBrsRWccS0SaQMAAAAAAAAA0EB4aQMAAAAAAAAA0EDOtuT3kYjY0xc/5MNQedP7778/ED9hIBnfi9casDkE5rEFMIftgHnsfJjDdsA8dj7MYTtgHjsf5rAdnHEezyqnDQAAAAAAAAAA9A/IowAAAAAAAAAAGggvbQAAAAAAAAAAGggvbQAAAAAAAAAAGggvbQAAAAAAAAAAGggvbQAAAAAAAAAAGggvbQAAAAAAAAAAGggvbQAAAAAAAAAAGggvbQAAAAAAAAAAGggvbQAAAAAAAAAAGsj/AzYwWoI30EbfAAAAAElFTkSuQmCC\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: 3.35 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n = 10\n",
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(x_test[i].reshape(28, 28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "n = 10\n",
    "plt.figure(figsize=(20, 8))\n",
    "for i in range(n):\n",
    "    # What does encoded images looks like?\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(encoded_imgs[i].reshape(4*1,1)) # 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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0171, PSNR (dB) = 17.67\n",
      "MSE = 0.0558, PSNR (dB) = 12.54\n",
      "MSE = 0.0053, PSNR (dB) = 22.80\n",
      "MSE = 0.0319, PSNR (dB) = 14.96\n",
      "MSE = 0.0210, PSNR (dB) = 16.78\n",
      "MSE = 0.0067, PSNR (dB) = 21.72\n",
      "MSE = 0.0591, PSNR (dB) = 12.28\n",
      "MSE = 0.0462, PSNR (dB) = 13.36\n",
      "MSE = 0.0892, PSNR (dB) = 10.50\n",
      "MSE = 0.0258, PSNR (dB) = 15.88\n",
      "Wall time: 4 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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6877\n",
      "0.3571\n",
      "0.9114\n",
      "0.7184\n",
      "0.6458\n",
      "0.9014\n",
      "0.2335\n",
      "0.4002\n",
      "0.0990\n",
      "0.6421\n",
      "Wall time: 305 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/imamun93/Image-Similarities-using-SSIM/blob/master/blog5_ssim.ipynb\n",
    "import numpy as np\n",
    "from tensorflow.image import ssim\n",
    "from tensorflow import convert_to_tensor\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = np.asarray(ssim(convert_to_tensor(original), convert_to_tensor(contrast), max_val=1)) # <class 'tensorflow.python.framework.ops.EagerTensor'>\n",
    "    print(f'{a[0]:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9175\n",
      "0.5774\n",
      "0.9613\n",
      "0.8687\n",
      "0.8122\n",
      "0.9764\n",
      "0.6746\n",
      "0.7032\n",
      "nan\n",
      "0.8636\n",
      "Wall time: 900 ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\public\\fyp-code\\plaidml\\.venv\\lib\\site-packages\\ipykernel_launcher.py:145: RuntimeWarning: invalid value encountered in power\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(\"./Deep Autoencoders in Keras - 784-256-128-64-32-16-8-4-8-16-32-64-128-256-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
}