Simple Autoencoders in Keras - 784-16-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 = 16"
   ]
  },
  {
   "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": [
      "1595258355\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 27:23 - loss: 0.6920WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (7.399499). Check your callbacks.\n",
      "WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.716004). Check your callbacks.\n",
      "WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.358505). Check your callbacks.\n",
      "235/235 [==============================] - 53s 226ms/step - loss: 0.3097 - val_loss: 0.2264\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 32s 137ms/step - loss: 0.2008 - val_loss: 0.1798\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 31s 133ms/step - loss: 0.1711 - val_loss: 0.1619\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 29s 122ms/step - loss: 0.1578 - val_loss: 0.1518\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 30s 127ms/step - loss: 0.1498 - val_loss: 0.1458\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 32s 134ms/step - loss: 0.1448 - val_loss: 0.1418\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 26s 113ms/step - loss: 0.1322 - val_loss: 0.1304\n",
      "Epoch 14/50\n",
      "235/235 [==============================] - 29s 124ms/step - loss: 0.1315 - val_loss: 0.1298\n",
      "Epoch 15/50\n",
      "235/235 [==============================] - 29s 124ms/step - loss: 0.1310 - val_loss: 0.1294\n",
      "Epoch 16/50\n",
      "235/235 [==============================] - 27s 115ms/step - loss: 0.1305 - val_loss: 0.1289\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 29s 124ms/step - loss: 0.1302 - val_loss: 0.1286\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 28s 117ms/step - loss: 0.1299 - val_loss: 0.1283\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 27s 114ms/step - loss: 0.1297 - val_loss: 0.1280\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 26s 110ms/step - loss: 0.1295 - val_loss: 0.1278\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 29s 123ms/step - loss: 0.1293 - val_loss: 0.1277\n",
      "Epoch 22/50\n",
      "235/235 [==============================] - 29s 122ms/step - loss: 0.1292 - val_loss: 0.1276\n",
      "Epoch 23/50\n",
      "235/235 [==============================] - 23s 97ms/step - loss: 0.1291 - val_loss: 0.1274\n",
      "Epoch 24/50\n",
      "235/235 [==============================] - 24s 103ms/step - loss: 0.1290 - val_loss: 0.1273\n",
      "Epoch 25/50\n",
      "235/235 [==============================] - 27s 116ms/step - loss: 0.1288 - val_loss: 0.1272\n",
      "Epoch 26/50\n",
      "235/235 [==============================] - 24s 102ms/step - loss: 0.1287 - val_loss: 0.1272\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 25s 107ms/step - loss: 0.1287 - val_loss: 0.1270\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 24s 101ms/step - loss: 0.1286 - val_loss: 0.1269\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 21s 91ms/step - loss: 0.1285 - val_loss: 0.1270\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 22s 95ms/step - loss: 0.1284 - val_loss: 0.1268\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 21s 90ms/step - loss: 0.1284 - val_loss: 0.1267\n",
      "Epoch 32/50\n",
      "235/235 [==============================] - 22s 93ms/step - loss: 0.1283 - val_loss: 0.1267\n",
      "Epoch 33/50\n",
      "235/235 [==============================] - 22s 94ms/step - loss: 0.1282 - val_loss: 0.1267\n",
      "Epoch 34/50\n",
      "235/235 [==============================] - 21s 90ms/step - loss: 0.1282 - val_loss: 0.1265\n",
      "Epoch 35/50\n",
      "235/235 [==============================] - 22s 95ms/step - loss: 0.1281 - val_loss: 0.1266\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 22s 95ms/step - loss: 0.1281 - val_loss: 0.1264\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 20s 84ms/step - loss: 0.1280 - val_loss: 0.1265\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 20s 86ms/step - loss: 0.1280 - val_loss: 0.1264\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 21s 88ms/step - loss: 0.1279 - val_loss: 0.1264\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 18s 78ms/step - loss: 0.1279 - val_loss: 0.1262\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 20s 86ms/step - loss: 0.1278 - val_loss: 0.1263\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 18s 76ms/step - loss: 0.1278 - val_loss: 0.1262\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 17s 73ms/step - loss: 0.1278 - val_loss: 0.1262\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 22s 94ms/step - loss: 0.1277 - val_loss: 0.1261\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 20s 85ms/step - loss: 0.1277 - val_loss: 0.1261\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 18s 77ms/step - loss: 0.1276 - val_loss: 0.1261\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 16s 68ms/step - loss: 0.1276 - val_loss: 0.1260\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 16s 67ms/step - loss: 0.1276 - val_loss: 0.1260\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 15s 66ms/step - loss: 0.1275 - val_loss: 0.1260\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 14s 60ms/step - loss: 0.1275 - val_loss: 0.1259\n",
      "1595258355\n",
      "Wall time: 21min 12s\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": 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+WH4yJAAATEElEQVR4nO3b24uVdf/G8XtyTMuaKcvUqJRKIirbYLQ7C6LwoCAiiEI8CKKTCLWDNgdBoBBU0F4izIQo2x+oQURQ0VOgUYhKIZSStBn3WFOps35/wK8H8l7fub/X9Lxeh7L6cNHqbs28WQ70er0GAAAAgCzH1R4AAAAAwP8n2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAoMFjefGkSZN6g4PH9I/EOOWUU2pPaO3gwYPN6OjoQIlbQ0NDvZkzZ5Y41bn9+/fXntCX3bt37+71ejNK3BoeHp6w7+OUKVNqT2ht165dzd69e4s8i6effnpv7ty5JU51bmRkpPaEvuzcubPYszg4ONibPHlyiVOd++OPP2pP6Euv1yvyLJ588sm9GTOK/OfQualTp9ae0Jdt27YVexZPPPHE3vDwcIlTnev1erUntHbw4MHm999//5//XNyxY0ftCX0p/TPqrFmzSpzq3I8//lh7Qmt//vlnc+TIkSLP4kR+D/ft21d7Ql9GRkb+9lk8pgIzODjYnHnmmeVWdejmm2+uPaG1119/vditmTNnNk8++WSxe116//33a0/oy8svv1zsE33mzJnN888/X+pcp84777zaE1or+f+RuXPnNhs3bix2r0svvvhi7Ql9uffee4s9i5MnT27OPffcUuc6tXXr1toTIsyYMaNZvnx57RmtXHDBBbUn9OWKK64o9iwODw83ixcvLnWuU2NjY7UntLZ69epityby5+Ldd99de0JfSv6MOmvWrGblypWlznVq2bJltSe0tm3btmK3Zs2aNWF/1nvzzTdrT+jLCy+88LfPor8eBQAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAg0ey4tPOOGEZv78+eO1ZVw9/fTTtSdEGB0dbbZt21Z7RitHjhypPSHG77//3mzatKn2jFa2bNlSe0JrBw4cKHbr559/bh5//PFi97q0Zs2a2hNizJgxo7nnnntqz2jl+++/rz2htddee63YrbGxsebQoUPF7nXp448/rj0hxuTJk5vZs2fXntHK5s2ba09o7ejRo8Vubdq0qRkYGCh2r0vvv/9+7Ql9efnll4vdOu6445oTTjih2L0uPffcc7UntLZ48eJit3bu3Nncd999xe516fDhw7UnjAvftAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAINDgsbx4dHS0+eabb8Zry7i6+OKLa09obfv27cVujY2NNYcOHSp2r0uvvPJK7Ql9Wb16dbFbkydPbmbPnl3sXpcWLVpUe0Jra9asKXZrz549zapVq4rd69LXX39de0Jfpk6dWuzW/v37m/fee6/YvS799NNPtSe0dvDgwWK3RkdHm61btxa716UlS5bUntCXpUuXFru1b9++5u233y52r0uffPJJ7QkRpk2bNmF/Xn/xxRdrT4jx559/Nt9//33tGa1M5M/Fkr/f9Xq95q+//ip2r0tjY2O1J4wL37QBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQKDBY3nxnDlzmueee268toyrrVu31p7Q2hNPPFHs1t69e5s33nij2L0uffTRR7UnxBgZGWmef/752jNamcjP4q5du4rdmjt3brNq1api97q0aNGi2hNizJ49u3nkkUdqz2jltttuqz2htcOHDxe79euvvzZPPfVUsXtd+uGHH2pPiDE0NNTceOONtWe0MjQ0VHtCa59++mmxW0ePHm1+++23Yve6tHDhwtoT+rJhw4Zit6ZMmdLMmTOn2L0u3XHHHbUnRDj11FObW2+9tfaMVvbs2VN7Ql+2b9/+t3/umzYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACDR4LC/++eefmxUrVozXlnH1xBNP1J7Q2sqVK4vdGhsba0ZHR4vd69JNN91Ue0JfPv/889oTInzwwQe1J7R24MCBYrcmTZrUDA8PF7vXpWeffbb2hL6sXbu22K0DBw4069evL3avS/Pnz689obVNmzYVu3X88cc3Z511VrF7Xdq8eXPtCTH27NnTvPrqq7VntHL99dfXntDa8ccfX+zW8PBwc+ONNxa716WJ/HtGadOmTWuuueaa2jNamci/a5T8PWN0dHTCfr7ceeedtSf05aWXXvrbP/dNGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAg0ey4vPOOOM5v777x+vLePqySefrD2htV9++aXYrVmzZjXLli0rdq9LK1asqD0hxrnnntu89dZbtWe0sm7dutoTWlu+fHmxW4cOHWq++OKLYve6tGXLltoTYoyMjDQrV66sPaOVZ555pvaE1nbs2FHs1kUXXdR88sknxe51acmSJbUn9GX79u3Fbk2fPr25/fbbi93r0pEjR2pPaG3SpEnFbp144onNZZddVuxel2655ZbaE/qyZs2aYrf279/fvPvuu8XudenKK6+sPaG1b775ptitmTNnNg888ECxe1368ssva08YF75pAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAA71e7x+/ePr06b0bbrhhHOeMnzfeeKP2hNYWLFjQbNy4caDEraGhod7VV19d4lTnPvzww9oT+rWp1+stKHFoeHi4d+2115Y41bl33nmn9oTWrrvuuuarr74q8iwODAz88//5hrn99ttrT+jL2rVriz6L1113XYlTnZs2bVrtCa199NFHzd69e4s8iyeddFLvkksuKXGqc//5z39qT+jLwMBAsWfxpJNO6l166aUlTnXu/PPPrz2htXXr1jW7d+8u8ixOmzatd+GFF5Y41bmFCxfWntCXxx57rNizODQ01LvqqqtKnOrcww8/XHtCa/fcc0/z7bff/s//jHosbSPRf/tc9E0bAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACDR7Li0877bTmrrvuGq8t42rPnj21J7R25MiRYreOHj06Yf9dPProo7Un9KXk/oMHDzYffPBBsXtd+uyzz2pPaO3QoUPFbp199tnN0qVLi93r0tjYWO0JfVm7dm2xW/PmzWvWr19f7F6XrrnmmtoTWvvrr7+K3frjjz+a7777rti9Lr311lu1J8Q455xzmmeffbb2jFYuv/zy2hNaW7BgQbFb8+bNazZs2FDsXpcefPDB2hNiTJ06tbngggtqz2hl3759tSe0dvTo0WK3zjnnnOahhx4qdq9LK1eurD1hXPimDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACiTYAAAAAgUQbAAAAgECiDQAAAEAg0QYAAAAgkGgDAAAAEEi0AQAAAAgk2gAAAAAEEm0AAAAAAok2AAAAAIFEGwAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQKINAAAAQCDRBgAAACCQaAMAAAAQSLQBAAAACCTaAAAAAAQSbQAAAAACDfR6vX/+4oGBkaZpdozfHP6LOb1eb0aJQ97DqryPE5/38N/B+zjxeQ//HbyPE5/38N/B+zjxeQ//Hf72fTymaAMAAABAN/z1KAAAAIBAog0AAABAINEGAAAAIJBoAwAAABBItAEAAAAIJNoAAAAABBJtAAAAAAKJNgAAAACBRBsAAACAQP8HUrGXafaiN9wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dW8xV1fW+B9Vai0A9AQVRxKqAiuL5DFVajWg0tTVGY2x70yZNY5p416YxNV6QNm2aHpJeWL1obGyrRtqmioZaD4hAEeUMghwU5eiJWtQe+F/8/p2+8+Vbk/3tb3+w9trPczU+5mTttedcY865V8Y7xpC9e/cGAAAAAAAAAADUi08c7BsAAAAAAAAAAIB94aUNAAAAAAAAAEAN4aUNAAAAAAAAAEAN4aUNAAAAAAAAAEAN4aUNAAAAAAAAAEAN4aUNAAAAAAAAAEANObQ/nYcMGUJ98IPE3r17h3TiOszhQWXn3r17R3biQszjwQNfbAT4YgPAFxsBvtgA8MVGgC82AHyxEfTpi0TaABw4Nh3sGwCAiMAXAeoCvghQD/BFgHrQpy/y0gYAAAAAAAAAoIbw0gYAAAAAAAAAoIbw0gYAAAAAAAAAoIbw0gYAAAAAAAAAoIb0q3oUwEAYMuTjhOaHHHJIsj/xifzdof596KGHVvbT6/33v//N2vbu/Tjpuf4//feIiI8++ijZ//rXvyqv4f8PAAAAAAAAYLAh0gYAAAAAAAAAoIbw0gYAAAAAAAAAoIYgj4J+o9Km0aNHZ22TJ09O9oknnpi1zZgxI9nHHXdcsj/96U9XftaePXv6tCMidu/eney33nora9O/161bl+w1a9Zk/V5//fVk79q1K2v78MMPk/2f//yn8h5VmoWMamCo5C1iX0nc/yjJ4QB6HfejKvAbgM6j+1are5qfMXyPA4D60epeq7DvQrsQaQMAAAAAAAAAUEN4aQMAAAAAAAAAUEN4aQMAAAAAAAAAUEPIaQMt8clPfjLZEyZMSPYVV1yR9bv66quTfdppp2Vtmv9G89hojpyIao2o60C1RPc777yTtWnums2bNyd7+/btWb8dO3YkW3PYROQaczSo/UfnUe1PfepTWT99LiZNmpS1HXbYYclevXp1sjUXUUSe74i5ah/Pt3DEEUckW+dp2LBhldfw/FI7d+5M9gcffJBscjYMHJ0vXUcPPTTf2tX/dN3897//nfVrx3fa0fS3+1mwL+6zumYeeeSRWdvhhx+e7LfffjvZ//jHP7J+pRxu8H/4c68+p+cbXUMj8rOUjrvm6IvI/WMwfKXVcxZAUyjll6raSyMihg4dmuyxY8dmbccff3yfn+W/NTZt2pTs9957L2vTfXiw/R66GyJtAAAAAAAAAABqCC9tAAAAAAAAAABqCPIo6BMPI9QwPQ0j1FDfiIgRI0ZUtmkIoIbol2QSpZBFvb6HIKsER6VTGhIeEfHRRx8l20PCCU3cPyVphLbp3HmJdy0NP2XKlKxN50fDTbdt25b1a7U8O+yL+thnPvOZrG3q1KnJ/tKXvpTs008/Peun/rx06dKs7ZFHHkn28uXLk/3Pf/4z64dcqm90flz2pL6kc6frcEQ+P+pHLotxuVTVfejaq3KciNzX9XP9b/0sn/teXXtblZppPx9/XU8vvfTSrE37Ll68ONkrV67M+r3//vvJ7tW5iNh3PvTZdonouHHjkq372HHHHZf107GdN29esteuXZv1K/li1T2W7tfXjqqy5O6zugf36rPQrgS0LvTqvEWU90/1D/3N4BKoG264Idm33HJL1jZq1Khkq++88cYbWb8nnngi2Y899ljWtm7dumSrdMrPsrpPcs7tTYi0AQAAAAAAAACoIby0AQAAAAAAAACoIQdVHtWqtKLU18P+2gkD7OXQwSpK46pVDjQjekTEwoULk/3mm29mbRpSv2vXrmRv2LAh66fVZrT6hUtnzjrrrMr7X7BgQbJXrVqVbA1NjqBCVCt0whc1rFPDrR3N0u+fpaGsnfD7XkbHVSU2Lnv66le/muxp06Yl22VUGnbsVeNOOOGEZP/qV79K9qJFi7J+Kpfq5fksVbgoVV5TWYzLRVUSpdW9fJxLUguV1hx99NHJHjly5L5f4v/jclT9W9eEpkvjSuOqf7daOaRKrhwRcfbZZyf75ptvztp0DdVnRCssRvS2L+p4arWtiFwKceGFF2ZtX/jCF5Kta55LEJ955plk61nH98XSuOszo+cql2zp/ftzomcftV22qjKtpkkySmttSX6vkje1XdKmEpktW7ZkbTrOpXHV56C0dqjt62nT5m1/6FiUUizonqaS4okTJ2b9vvKVryT75JNPrryGjrtLlHXvdpnW3/72t2Tr7xWvxKnyq16b04FS2oN1PnTN9POWrtEuI61aTzuddoNIGwAAAAAAAACAGsJLGwAAAAAAAACAGsJLGwAAAAAAAACAGjIoOW1KZQdV/+f6W81nUSpPqLpB14dp6V+9nusLVaOq/yeiugy061U1L0uplHQTNOH6HfS7el4KLR3qGm3V8Krt2kD9LH12VCceEXHjjTcmW/M5RES8+OKLydZ8DqXPgo8p6T/boZSzQvMW+RyrJlzzETGPA0P9SvPY3HXXXVk/Lfk9fPjwZPszoeOv8xkRce655yb7tttuS7avp0uWLEm251Xo5fnVvfCoo47K2jSvl+YS8vF77bXXkq1+VNLFl3LaaB4HzzWm11y2bFnWpvnQemlOq3IsOK3msanKSRURMXPmzGSfd955WZvOjeaj0/wdvYiOp/qb5m6KiLjsssuS/cUvfjFrmzx5crL1TLlmzZqs39y5c5O9Y8eOZPcnr5M+Q7ouuy/qeqHnoIiIbdu2JVv90p/PqrN3RHf6cFWevIj8N8nnPve5ZPtcX3rppckeP358snUcIyJeeumlZD/00ENZm55RdU1u9zlQur1E+UDRcSn9DqzyI/89oTn8fGx1vtTv33nnnayf+pv/XlQ/0rOtn3P97yZRem+ge5zPobbpuWTMmDFZP/1t4X6j86bXP/bYY7N+emb1fK561lHb8xL5ube/EGkDAAAAAAAAAFBDeGkDAAAAAAAAAFBDBkUe5dImDV/S0Pnjjz8+66flFD0sqap0oYf16jVPOeWUZI8bNy7rp2FmGzduzNq0RKN+roc1aXjUH//4x6xNQ/1dftWNaPiefh8P/dK/PfS+JGtrBS+/edJJJyXbw5g1lF/vtwnhvd2Ihj66LFJL1OqcRuQh3Rpuijyqf3joqYaOfve73032RRddlPXTtbtUUlTnw2WRVVKsa6+9Nuun667KLPu6ZpPx0F2V8qrULCLiuuuuS7aGd8+fPz/rp+tyq9Jdb9N5VEmA7rMRuTRr7dq1WdsHH3yQbN0fmua/pdK8TqvfvUpipXMREXHJJZckW8P6I3I/VR9+77332rqnJqLfXc8REbmP+RlV/9+bb76Z7Keffjrrp+fGVkPl/fnRUrRnnHFGsqdNm5b107P4ggULsjYt817ll319drej38dL+qq8QiVQ6lMREaeeemqy9RnZs2dP1k/HX2U6EftKs6oo+WKVvLXX/Nd/c+rvPZW5+W8x3RfVB1SWHxGxYcOGZPszo2uxymL+9Kc/Zf30TKNy5YhcJlmSLzdtXnXsXOar8zZjxoxkX3nllVk/lYTq+lySQG3ZsiVr27p1a7L1vOVydJ37d999N2vT+3/llVcq72OgEGkDAAAAAAAAAFBDeGkDAAAAAAAAAFBDeGkDAAAAAAAAAFBDBiWnjaPaTdV1uoZNS3K5lkz116ol836aM0f1ba5D1HLdrl/U3Cknn3xysj1vipZY9fJuK1asSHaTy3+7Jlvb2v2uqjlWjfGPf/zjrJ/mQnF94VNPPZVs1ar2p5wiDIyqMqqf/exns35aVtpzMGhOjJ07dya7VKoY9sX19N/+9reTrWVsvZ/Oofqz+lREvob6eqq+qdcbO3Zs1k/LqupcR+Qa5G5aQ6vGz1Hds+9VZ555ZrJvvfXWrG3SpEnJfvXVV/u0I3LNfLu+o/njNP+c56bTPAGeQ0D3i26axwNFq2Oic6Fl3yPyXCuup9d8CbNnz0625pPqz300harv6+cF3Z8094H31Xwxmt8gIl8725nviPzs8/Wvfz3ZxxxzTNZP82ho3oaI/Myqfur3NNBchAcbz8lTVeY5Is+pp+ua51PUXHs613//+9+zfroX+hlVz0St7hPe1o3z0Sl0/PQMGVF9lvC8Tq+//nqydV/861//mvXTHFWaQyUiYv369cnWMu4+3yU/0rYm/0ZxX9R3AOeff37WdttttyVbc/l5OXY9L+lvbc/TtmbNmmS7n+q+qJ/lv1X0fYPn1dR3CrqfDrTEt0OkDQAAAAAAAABADeGlDQAAAAAAAABADRkUeZSHX2s4qIaMeXiplkJcvnx51qYhbno9D7dSCZOWDPOScKtWrUq2S5s05Fyv4VIsDY8aPXp01uaf1yT6E77ZDlrCdO7cucn2sLgqyVxExOrVq5PdTjhyRP5slcJXezlEtdXvrmUxJ0+enLXpvGqocUTEokWLkq2lhHt5zFtFZalelvmqq65K9ogRI5Lt62mVJGrevHlZv/vuuy/ZLo/SENgLLrgg2VqqNiIvV+uhrb/97W+TvXv37ugWWi2Xq3PlMtzLL7882V7eWUNvNUx7yZIlWT8tRduu7+j/Uxmy730aouwyuiaHfrdKu/ITlXWoXMPDylVKo3MRkcsDtBQ1ctOPUZ91WcyYMWOS7TJGnTsNt9ewee9XKgere6aeSSMi7rzzzmSfc845yV63bl3Wb+nSpcnW83VEtWy86ecbHVeXk+l5X8tGeynvxYsXJ/uZZ55Jtspo/PqeDkL3Xd3vSmWemzYX/cFLpOv54fbbb8/ahg0bluwXXngh2Sq3j9hXvvs/9KwZETF//vxku8RKfwe2ur+1ei5oGv67WM8zWtY7IpeA6rrrc6byU93TXOKmvy10fY7IZVAqi3Tpsa7XfkbVZ2swU6IQaQMAAAAAAAAAUEN4aQMAAAAAAAAAUEMGRcPj4UAawq1ZlT2bvUqgPDxNQ+JL2Zg17ExDQz0sqyosLiKvoKLVAryyisqqXBLQtCoZg/kdNEw0IuKBBx5ItkpnPJRY59DlGpodvNXQ71KVAaUJ8zkQ2vn+6n8TJ07M2tRXPKu7hjQiregfKpPw0FOtyqZhxz63Kjt89NFHk/2Tn/wk66dVGPwaGqqs1fg0FD0iD4FV+U1ELpN7+eWXk113WUervqJrjcuGdR5dbqShwY888kiyfT9qZ5xKMh6VCruca+PGjcn2MORe9eFO7xkq23HJnM61S3N+/vOfJ7ubZIaDje796ou+Dqkv+nlQ51j7uayjSjrvFaJ0jb7pppuytosvvrjPe1fZf0Qu5fBw/qo1oWnyKD/XqTxKK9ZG5BVhVVK8bdu2rJ9W5VJJlEvmVHbhz5I+F/pbwqX+na4+0634vqhnGq2wGBHx7LPPJlsrOvn+2SpeKWiguE+1moqhG9Hvo74XkfufSucj8iqIKjd6/vnns36zZs1K9oYNG5Ltc6Zj6X6qUkhdA9xn9ZoutdNzz2D6LJE2AAAAAAAAAAA1hJc2AAAAAAAAAAA1hJc2AAAAAAAAAAA15IDktFEdu+poNb+N/z/XHqqWrFWdn35WqcSa6+y0hKbmWXA9vuqHVTfp9wv7orruW265JWs76aSTkq1z6OVLtYTlvffem7V5TqRW8Bw2eo91z51RR9TnNDeUzm9E7isrVqzI2trVIPcivsYNHz482ZoDISJ/tjU3lOamiYi46667kj1nzpxku+6+VCZXS0Jff/31yfacNroOay6riIhLLrkk2ZrTpu6U9qoqHbvPo/qHl+194oknkq2+42tlO9p4vw/Vd2u+B8/Robry0h7fZNr9nqVziu5PqsEfO3Zs1k/3Ktfda86wXpmL/qK5ZbxMcwnNUaL7nc+PzqPmVtD8JxH5OVTXTb/++vXrk/3cc89l/bZv357s0hmmyc+C+5TuT55fb+rUqcnWfBv+G0H3Iz3PaB6OiLxUu+fH0L1W11DPz6i5iPw3SJPnLSKfu1GjRmVt6h/up8uWLUu251WrO6WzQDfmhCudDdX/tNR2RL5O6jr20EMPZf0072Upl4x+tudAuvvuu5OtZb59/Lds2ZLsl156KWvTc/Rg5rQl0gYAAAAAAAAAoIbw0gYAAAAAAAAAoIYMijyqhIZ3eQi3hhR1Qo5SCkvSEFgv+3fttdcmW0NeNTQqIuLxxx9PtpZejUBO0xcaaqYhw1oGOCIPG9XQUJc8aZicy2qq5r5U1ttDYHXumc+P0TErhfNryK+GPmr4YUQ+tjt37szaStdvFb1ffS6aFlrsY6VlmV3CsmPHjmSrhOIHP/hB1k+lSCXJp/qKl1M899xzk62SKA9p1mt4m64ROp/d7Jc6X/rdPYRY18Dly5dnbTp3vp8OFL+PK664Itnqzy6B0r2wHZlqr1EVDu9yXfUJlWdoaeiIfF1buHBh1sZ89I2Oe1VJ7oi8rKvLnlS+oXN1+eWXZ/30+ipzGzlyZNZP5VLaLyL3dQ3T95LfpTW76lnzfbGb19iI8pnPy0gfffTRydY59PFXuZpef9iwYVk/3Qvdn3UvVLmxS5R1jXe5eNPOMI76ipZijsh/t/nYDh06NNn6W8+f5bqMX0kS1e3o9/FznUr4/fyia9fu3buT7evkqaeemmx9h6C+HBFxzTXXJPvGG2/M2lTCr/frZ5s33ngj2W+++WZUMZjPFZE2AAAAAAAAAAA1hJc2AAAAAAAAAAA15IDIozRUqJTdWaVTrVbdKH2W4iGvGnp63XXXZW0adqdhWZ6Z/9lnn032O++8k7V1e0hpJ/B50tA4DWnbunVr1k9lZxrm6MydOzfZe/bsqeynoZMeRqnhq/5ZOoetVjFSmYNfo1vxedTvWAqr1vG87LLLkj1hwoSsn0o+1N9K91FaH3yOmyKn2R/+7Oka56GnKpt44IEHku0Z8UvrdRVe4UB9U9v6I1XsTxWXbqFKnuDPr/qEP7+6r6m/uSyiqpqBz4GG91966aVZ22233ZZsfba8mpeG9/vzU5dw9AONfm+fX/VbnU/3Z5XjqIR7xIgRWb+33nor2StXrszaOrH+tboOdys67ipNjIhYvXp1sl0+WCXXcFmHVhjyikKKyq38Wdi2bVuyFy9enGw/S+l6689dlSTDn5Fun2+/Z10bfbxUnq3yCpXYROTnRt2rStI6R6+pclNP16D36HOj36UbKwvtDx3n8847L2tTmYzPz/Tp05Ot8rI1a9Zk/XT89BpaNdh57bXXsr/1fNNqlWPfd1uVZ3b7+dW/m46dj6tWyNNx9fVUZYb6TGjltoi84qXLGHVt1M/SNAIRuTzWf6voGXswfZFIGwAAAAAAAACAGsJLGwAAAAAAAACAGsJLGwAAAAAAAACAGnJActpU6btKmj9va7UMWlUZVdcOT548Odlaktb/n+bbeOSRR7J+69evT7bnPGmivrS/uH5RtbqqL3SttZZVU9245w3Skmul8S6VnFMNpOe00XJvWiq19FndqPnuFD6Pqgm/+uqrk+1zoD7mZfSqcnGUPtufuybPSUkLrTkRvHSh6oc1J0knNNMffvhh9rfOr+Y7UR1/RD5P7mPt5NapG6VcSzqPXpZZtdOev0R127pWetlY1Y7rc6Jl4SMiTjzxxGTfeuutWdvZZ5/d5zXef//9rJ/mVGEfbB/P63TGGWf0aftzpb69YsWKrK1UArqKUv6FpqytVefLt99+O/tbSzNv2bIla9P1TM8V7rN6vtT50LKzERFXXXVVsv18uWHDhmQvWrQo2e6LpfxVeh/az+e02+fY71/HSPNSRuT5hmbMmJFsXyf1rKj5MbyEuPqm72Gal2jTpk3J9t8qpdx0ujc0ca3VsfUcJepXfvaZNm1aso877rhk+28Iveb48eMrr6d72iuvvJK1Pfnkk8m+7777kq35kSLK+fyaVuZbUf/zPDBPP/10sj3fkK6HOh+6Bkfk/lG1Bkfsm/dIUb/SHFLz58/P+s2bNy/Z6rMR+bmNkt8AAAAAAAAAAD0GL20AAAAAAAAAAGrIAZFHKe2G8FWFG5XCzDS8WKUCERFXXHFFsk8//fSsTcOcZs+enWwvbaqSg24vxdYpNATtpJNOytpuv/32ZJ911lnJ9nA3DT/UkmubN2/O+nkosKKhvxqyquXhIvLQSUdlWlr6059h/buJz4H7XpUPewnUKVOmJFtLvHvI+VNPPZXsXbt2ZW1V41kqF+1t3R7e3SoupzjmmGP6tCOqS122I0P1v/35UF8shajqXHsYrZYu7tb5LMkO9Lt7CLeG4WoZzIiIiRMnJvuiiy5Ktpe01HHX9dZDfPX/qQQnImL48OHJVrmV36+2detcDSal50Btl1qoXEOlAS6d0RLQKuGOaP38VVUOutSvKXOtMhZ/tkvlotetW5dsXYtdvqbX1zaVH0bkZ1S/hkpr9Fzk8reSJKNKjto0mY0/lyrf9TPlb37zm2Q/99xzyR43blzWT2VUY8eOTbb/ztBzo3+WS53+R6mUvEtn1fdblZJ3E1ry231Ax8/bVL6mEijdwyLyddSvoaikX+c7Ii9Frp/7s5/9LOuna4n7WNVZoAnod1N/iIhYtmxZn3ZE9b6j58mI3Df1nOOlwdVn3Y/0N+fvfve7ZLu8+NVXX022p3LwtACDBZE2AAAAAAAAAAA1hJc2AAAAAAAAAAA1hJc2AAAAAAAAAAA15IDntOk0rZYN99wl06dPT7aXZHz++eeTrSUBtexbRHvlM5uGl8bTPDbf+ta3sraZM2cmW3W7rrvXv7WsppfY1PH3+1ANqpYX11LvEXm5OH+WVE+r+YtcU6naWs+z062likt+VaV593FRLbG2rVq1Kuun2vEqnbdTKp9copSfoRt14Pp99HmNyLW+rsPWvzWnl5epVB/TeffxVv/ztfamm25K9oQJE/q8XkSuCda8HBERc+bMqfx/3UIpN5SOs+fz0XXPcwJpDgXNGzVmzJisn65fCxcuTPaCBQuyfroeXnPNNVmbrmW6znkeKv0u3ehTg03pOVDb51rLsWveA837FhHx4IMPJtufpXZyA7a6ZnZzLjEdd12H/DupD7Sae6KUw0j3RV97tfSzr7cvvvhisjVXht9TKaeNtnXTXA0U/d5+3lBfUt/REusREUcffXSyNc+Y57V47bXXKj9LSxLruVRz2Hg/z6Ok96ufXcqZ0k3ovuXnRvUP/y3meb7+x5lnnpn9rTlo9Azj46X+576ovzWuvPLKZD/88MNZP32e3E+7dX5aQb9b6fdQq2Pg69jrr7+ebPU3Pw/rM7J8+fKsbdasWcnW3Dqe+0bPPe7PB2oOibQBAAAAAAAAAKghvLQBAAAAAAAAAKghXS+PKqHhUR7qraWfvezpX/7yl2SXwht7FQ0PPOqoo7I2lUJcddVVWZuGeVaVvfQ2DfN0CZSGR3pomsoDTj755GRriLm3+X1oGJ6W2PQQPw1L9bBHDe9sClVhgBpqGpGHourYuvRFQ3xbDTEshZyXaFqJ2pI8QcOs1fci8vlQP/Uyhhqyr+GlXpZUy9PqGhARMXXq1GTrM6KloSMiXn755WTfc889WZuu0d0qj3Kqwob9+6n8wfcqlTfpWuOlTZcsWbLf/xORy+Zcjqprp4Z6b9++PeunIcVN8LFOU5IRqe3rqe53+rx4WdKXXnop2YPhK1XSqW6ea7133c9LZ4JWS7eX0OtNnDgxaxs9enSyXQKn/qxSkJJkq5vnp7+UJH1Kq5JVZ9iwYcnWMdfSwRFlWb3KavSZ88/V6zetHPT+ePfdd5M9e/bsrE3TVegeGRGxcePGZOs54/zzz8/63XDDDck+5phjkq2St4hc8u2yVUXXC3+2dO7aPb92O534nn4Nnas77rgj2S7TV2nTvffem7VpShT12TrKDIm0AQAAAAAAAACoIby0AQAAAAAAAACoIY2TR2lY5GmnnZbsG2+8MeunUhuvoKEVozwbPOQhgBpSGBFx1llnJVsz7Efk4aE6rocddljWTyuYaAWi8ePHZ/00RN/D2DRkTitaqR2RywG88tPatWv7vHenJG1oIupjGio6ZcqUrN+FF16YbJXTrF69OuvXThW2Ulizh7S3E7beLej38Uz3GjLs4eI6H5///OeTrSH5fo1jjz022V6FQUNR3Z8VlZiuXLkya/vRj36UbM3gH9HMSn1Vz6U/oyp/0HDxiLyqydtvv51slz2pzFev55/13nvvVX6Whpmrv3mljl4L4R8oul6p7+heGpHvizr+8+fPz/r5OtAKpXWxm6tCtUNpD29VdlNCr6HnlK997WtZP5X3v/HGG1mbVkzB3/pHaQ6rzjZefVGrIKqs18+QpXOJtukcuv/q+urXr0rZ0BQf1bO1nxtL31Fl3TqeLqPSilT6W+a6667L+l1//fXJ9nnUs4nOj8u/e7VaW6fxqlC//OUvk33BBRck2/1cn58nnngia9NzT93nhkgbAAAAAAAAAIAawksbAAAAAAAAAIAawksbAAAAAAAAAIAa0ricNpqrQTXCWi4zImLr1q3Jvv/++7M21/JDjmoFvfy1luHzfEBVuS68bPg555yTbC3R5yW/NaeN51VQXaKWVvR70HvcvHlz1qYaSC2F7Jpj1bH6eDQB14Zqfh/Nc3LzzTdn/apKBHsp4U6g810q09fk/DauoX700UeTffHFF2dtOjdqe84n1Q+r//WnFK76hOZmuPvuu7N+zzzzTLJ7OZdYqVSor3O63qqPqe3/r1QmXue7lINBcyn4etiJvB9NprSejhs3Ltlf/vKXs356htH5Xbp0adav0+ta6XpNXk8jBuc76To6Y8aMZHs5Yu3nPqbrYxPHvZO0msPG/1bbc4RpPrZNmzYluz/5kNTvtc1zFOnvEc29EZHnU2l6PkX3gTVr1iT7iCOOyNqq8kx6ThvtV8qbp/udl/zW+9JnwedK6bU8YQNFfcV/Z8ycObPPfj6fP/zhD5O9ffv2rK2bxp9IGwAAAAAAAACAGsJLGwAAAAAAAACAGtL18igPVZs2bVqyNfTUQ6X+8BAv+b4AAAfiSURBVIc/JFtL9kV0V6jUgaIqbNQlGS+88EKyTznllKxNyyZqGL5K2ryfSptKpfb8PjQ0UUMWX3311ayfls5cvnx51rZixYpka3ish2lqiGUTy2+WykVPnjw52Vpuz//ftm3bkr1r167i9Vuh1ZD9/rR1Oy7NW7BgQbIffvjhrO2b3/xmsseMGZNslw+6z7WCh2lruWn93Hnz5mX9mljWuxPoM+tz7GH7Vf2qnnsNJ47Iy0qPGDEia9O1TX3YQ85LMromro/9pbSenn322ck+44wzsn46ljrmWup9MOjV9TSif6XQq/Yx97Hjjz8+2d/4xjeSreW/I3JJhpf8dplkO1Tdb9PmtD/fR/cuHeONGzdm/XQda1USX5KiapvPrfq6n3P1Ppo2b46fKzQtgcupda1U2/cf/f2o56BRo0Zl/UoSZU2d8Oc//znZLlGGMiX5oO6F99xzT9ZP51B94Pnnn8/66Rm4m32FSBsAAAAAAAAAgBrCSxsAAAAAAAAAgBrCSxsAAAAAAAAAgBrSNTltVO+meRemTp2a9fve976XbNUlau6SiIhf//rXySaXwv5RDaDqO72s3dNPP51s1XpGRIwfPz7Zmu/Gc99MmDChz3tQDXBEru/1MtLz589P9mOPPZbs9evXZ/20nKJrVavKKTa9tGJEdQ6jiNz/tFy7lh+OyHXgc+fOTbbr8zsxnt2sUR0s1D9+8YtfZG1aJnjWrFnJnjRpUtZP9cKl3EPqO+pvERHf+c53kq3rcC/4UacplbNv5d8j8nn0nHCa08bzBGzdujXZmhvM90/VopfKhveqz/qYaLnaE044Idmef0FzqWmeqKq8RgOhVBa+V+ctonW/Uh/wXDXXX399sjVXgz8XmhNj7dq1WZvmUdHP7eW5UUrPb1U//7u0P7U6zlXPRESer7FU3l3PqK3mz+kFSrneFJ0D36u0TXPk+Pxq3rCdO3dmbb///e+TrXlUNCdVBHvf/vD179hjj032nXfe2ee/R+TjqnNz6623Zv2a4jtE2gAAAAAAAAAA1BBe2gAAAAAAAAAA1JDayqM8VGrYsGHJnj59erK///3vZ/3OPPPMZGt4sZeXdekOtI6Oq4dma2jn5s2bs7YqyU0pbHTkyJHJ1pJ8EXm42/bt27M2DeXXe/KQc8IU+0bHxcdIQ0xVbvbggw9m/TTkV/3Pw0uZg8FBx1VDfyMi5syZk+xnn3022aeffnrWb+bMmckePXp0shcvXpz1mz17drK9pDtlngePgfqOSwd07pYsWZK16Zq6atWqZLv0WGVVrcq5mk5JoqFnnR07diR75cqVWb+hQ4cm+6mnnkq2+1unx7hX52x/tFoK3fvps6DnFpdTqI/5s9CJUP9WJUTdTmn8+/P/2kF926Woeu7VM6qX9da9m/Nr37Tqi74f6e8XPQf5mqppAHSNjsilwlqevSlynMFE/UNlwhERN9xwQ7KnTZvW5/+JqE4DoOeVJkGkDQAAAAAAAABADeGlDQAAAAAAAABADamVPErDng4//PCs7dRTT032HXfckewpU6Zk/bTCkGZdX7ZsWdaPyiWdwcMSOyGF0HA3rUi0Zs2aAV8b+o/PqYaULl++PNk+P+pjKpkgxPfgo3Oj87lgwYKsn/8NzcJ9UeUaixYtytqOPPLIZKskysOQdf1mn/0/SpVOdMwff/zxZL/88stZP50rrR7lFRwZ84NDlSTKZU8qFdaqeyNGjMj6qbxcpVJ+zU7sn726B3e6MppfT/92v1QpjfYrnaFdGtKr89YpdKx1Pl544YWsn467jzlVodpHx9XTX1x44YXJHj58eLK9ApieRTRFQ1P3QSJtAAAAAAAAAABqCC9tAAAAAAAAAABqCC9tAAAAAAAAAABqSK1y2lTlwIjIddslrZrq3TZs2JBsz81AOTaA9lD/U02+2gBQT0r5NjRXyrZt27I2LVlbKieveQLQ+O+Lj4meWTQ/kI9/qYw01IuSj61duzbZ6m+ajzEiP6N6meFO57TpJQZzvPzaOoeaOy4iL/Ot83vYYYdl/arOWzB4+G/MpuZHOdjouHqumkMP/fj1hPqOr6c//elPk62/+ZsKkTYAAAAAAAAAADWElzYAAAAAAAAAADWkVvKoElryS/GQw0MOOSTZK1euTLaWT4zoTGlqAACAbqUk1fFwZRh8kEA1D/ej3bt3J1vPr14uWqU1PAuDw4EcV/8s/Q2yZ8+ePu2IvCwyMh1oEvo8qzQ4IuL+++9P9pNPPpnsTZs2Zf0WLlyY7F74XU+kDQAAAAAAAABADeGlDQAAAAAAAABADeGlDQAAAAAAAABADRnSH03nkCFDDpqwVnWdWhqxVJpN9W3drgXdu3fvkP332j8Hcw4hFu/du/e8TlyIeTx44IuNAF9sAPhiI8AXGwC+2AjwxQaALzaCPn2RSBsAAAAAAAAAgBrCSxsAAAAAAAAAgBrS35LfOyNi0357DQIqb/KSeD3A+A5e66DNITCPDYA5bAbMY/fDHDYD5rH7YQ6bAfPY/TCHzaDPeexXThsAAAAAAAAAADgwII8CAAAAAAAAAKghvLQBAAAAAAAAAKghvLQBAAAAAAAAAKghvLQBAAAAAAAAAKghvLQBAAAAAAAAAKghvLQBAAAAAAAAAKghvLQBAAAAAAAAAKghvLQBAAAAAAAAAKghvLQBAAAAAAAAAKgh/w/ajujlKWPCxwAAAABJRU5ErkJggg==\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: 12.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n = 10\n",
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(x_test[i].reshape(28, 28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "n = 10\n",
    "plt.figure(figsize=(20, 8))\n",
    "for i in range(n):\n",
    "    # What does encoded images looks like?\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(encoded_imgs[i].reshape(8*1,2)) # 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.0114, PSNR (dB) = 19.45\n",
      "MSE = 0.0383, PSNR (dB) = 14.17\n",
      "MSE = 0.0052, PSNR (dB) = 22.84\n",
      "MSE = 0.0161, PSNR (dB) = 17.94\n",
      "MSE = 0.0167, PSNR (dB) = 17.78\n",
      "MSE = 0.0040, PSNR (dB) = 23.99\n",
      "MSE = 0.0239, PSNR (dB) = 16.22\n",
      "MSE = 0.0379, PSNR (dB) = 14.21\n",
      "MSE = 0.0352, PSNR (dB) = 14.54\n",
      "MSE = 0.0163, PSNR (dB) = 17.88\n",
      "Wall time: 28 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": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7646\n",
      "0.5488\n",
      "0.7849\n",
      "0.8468\n",
      "0.6899\n",
      "0.8354\n",
      "0.6055\n",
      "0.4570\n",
      "0.6227\n",
      "0.8093\n",
      "Wall time: 3.64 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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9498\n",
      "0.8628\n",
      "0.9825\n",
      "0.9292\n",
      "0.9434\n",
      "0.9900\n",
      "0.8909\n",
      "0.7736\n",
      "0.8700\n",
      "0.9498\n",
      "Wall time: 2.22 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/tensorflow/models/blob/master/research/compression/image_encoder/msssim.py\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from scipy.ndimage.filters import convolve\n",
    "\n",
    "def _FSpecialGauss(size, sigma):\n",
    "  \"\"\"Function to mimic the 'fspecial' gaussian MATLAB function.\"\"\"\n",
    "  radius = size // 2\n",
    "  offset = 0.0\n",
    "  start, stop = -radius, radius + 1\n",
    "  if size % 2 == 0:\n",
    "    offset = 0.5\n",
    "    stop -= 1\n",
    "  x, y = np.mgrid[offset + start:stop, offset + start:stop]\n",
    "  assert len(x) == size\n",
    "  g = np.exp(-((x**2 + y**2)/(2.0 * sigma**2)))\n",
    "  return g / g.sum()\n",
    "\n",
    "\n",
    "def _SSIMForMultiScale(img1, img2, max_val=255, filter_size=11,\n",
    "                       filter_sigma=1.5, k1=0.01, k2=0.03):\n",
    "  \"\"\"Return the Structural Similarity Map between `img1` and `img2`.\n",
    "  This function attempts to match the functionality of ssim_index_new.m by\n",
    "  Zhou Wang: http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "  Returns:\n",
    "    Pair containing the mean SSIM and contrast sensitivity between `img1` and\n",
    "    `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  img1 = img1.astype(np.float64)\n",
    "  img2 = img2.astype(np.float64)\n",
    "  _, height, width, _ = img1.shape\n",
    "\n",
    "  # Filter size can't be larger than height or width of images.\n",
    "  size = min(filter_size, height, width)\n",
    "\n",
    "  # Scale down sigma if a smaller filter size is used.\n",
    "  sigma = size * filter_sigma / filter_size if filter_size else 0\n",
    "\n",
    "  if filter_size:\n",
    "    window = np.reshape(_FSpecialGauss(size, sigma), (1, size, size, 1))\n",
    "    mu1 = signal.fftconvolve(img1, window, mode='valid')\n",
    "    mu2 = signal.fftconvolve(img2, window, mode='valid')\n",
    "    sigma11 = signal.fftconvolve(img1 * img1, window, mode='valid')\n",
    "    sigma22 = signal.fftconvolve(img2 * img2, window, mode='valid')\n",
    "    sigma12 = signal.fftconvolve(img1 * img2, window, mode='valid')\n",
    "  else:\n",
    "    # Empty blur kernel so no need to convolve.\n",
    "    mu1, mu2 = img1, img2\n",
    "    sigma11 = img1 * img1\n",
    "    sigma22 = img2 * img2\n",
    "    sigma12 = img1 * img2\n",
    "\n",
    "  mu11 = mu1 * mu1\n",
    "  mu22 = mu2 * mu2\n",
    "  mu12 = mu1 * mu2\n",
    "  sigma11 -= mu11\n",
    "  sigma22 -= mu22\n",
    "  sigma12 -= mu12\n",
    "\n",
    "  # Calculate intermediate values used by both ssim and cs_map.\n",
    "  c1 = (k1 * max_val) ** 2\n",
    "  c2 = (k2 * max_val) ** 2\n",
    "  v1 = 2.0 * sigma12 + c2\n",
    "  v2 = sigma11 + sigma22 + c2\n",
    "  ssim = np.mean((((2.0 * mu12 + c1) * v1) / ((mu11 + mu22 + c1) * v2)))\n",
    "  cs = np.mean(v1 / v2)\n",
    "  return ssim, cs\n",
    "\n",
    "\n",
    "def MultiScaleSSIM(img1, img2, max_val=255, filter_size=11, filter_sigma=1.5,\n",
    "                   k1=0.01, k2=0.03, weights=None):\n",
    "  \"\"\"Return the MS-SSIM score between `img1` and `img2`.\n",
    "  This function implements Multi-Scale Structural Similarity (MS-SSIM) Image\n",
    "  Quality Assessment according to Zhou Wang's paper, \"Multi-scale structural\n",
    "  similarity for image quality assessment\" (2003).\n",
    "  Link: https://ece.uwaterloo.ca/~z70wang/publications/msssim.pdf\n",
    "  Author's MATLAB implementation:\n",
    "  http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "    weights: List of weights for each level; if none, use five levels and the\n",
    "      weights from the original paper.\n",
    "  Returns:\n",
    "    MS-SSIM score between `img1` and `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  # Note: default weights don't sum to 1.0 but do match the paper / matlab code.\n",
    "  weights = np.array(weights if weights else\n",
    "                     [0.0448, 0.2856, 0.3001, 0.2363, 0.1333])\n",
    "  levels = weights.size\n",
    "  downsample_filter = np.ones((1, 2, 2, 1)) / 4.0\n",
    "  im1, im2 = [x.astype(np.float64) for x in [img1, img2]]\n",
    "  mssim = np.array([])\n",
    "  mcs = np.array([])\n",
    "  for _ in range(levels):\n",
    "    ssim, cs = _SSIMForMultiScale(\n",
    "        im1, im2, max_val=max_val, filter_size=filter_size,\n",
    "        filter_sigma=filter_sigma, k1=k1, k2=k2)\n",
    "    mssim = np.append(mssim, ssim)\n",
    "    mcs = np.append(mcs, cs)\n",
    "    filtered = [convolve(im, downsample_filter, mode='reflect')\n",
    "                for im in [im1, im2]]\n",
    "    im1, im2 = [x[:, ::2, ::2, :] for x in filtered]\n",
    "  return (np.prod(mcs[0:levels-1] ** weights[0:levels-1]) *\n",
    "          (mssim[levels-1] ** weights[levels-1]))\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = MultiScaleSSIM(original, contrast, max_val=1.0)\n",
    "    print(f'{a:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Simple Autoencoders in Keras - 784-16-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
}