Simple Autoencoders in Keras - 784-4-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(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": [
      "1595258110\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 16:55 - loss: 0.6932WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (4.943012). Check your callbacks.\n",
      "WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (1.180006). Check your callbacks.\n",
      "WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.590505). Check your callbacks.\n",
      "235/235 [==============================] - 45s 191ms/step - loss: 0.3711 - val_loss: 0.2542\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 25s 105ms/step - loss: 0.2466 - val_loss: 0.2349\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 24s 103ms/step - loss: 0.2278 - val_loss: 0.2205\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 23s 97ms/step - loss: 0.2197 - val_loss: 0.2161\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 20s 85ms/step - loss: 0.2163 - val_loss: 0.2134\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 21s 91ms/step - loss: 0.2138 - val_loss: 0.2114\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 19s 80ms/step - loss: 0.2119 - val_loss: 0.2098\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 23s 98ms/step - loss: 0.2103 - val_loss: 0.2085\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 19s 80ms/step - loss: 0.2091 - val_loss: 0.2074\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 27s 116ms/step - loss: 0.2019 - val_loss: 0.2007\n",
      "Epoch 26/50\n",
      "235/235 [==============================] - 27s 114ms/step - loss: 0.2018 - val_loss: 0.2005\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 27s 117ms/step - loss: 0.2016 - val_loss: 0.2004\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 26s 113ms/step - loss: 0.2015 - val_loss: 0.2003\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 26s 112ms/step - loss: 0.2014 - val_loss: 0.2001\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 27s 114ms/step - loss: 0.2013 - val_loss: 0.2002\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 25s 107ms/step - loss: 0.2012 - val_loss: 0.2001\n",
      "Epoch 32/50\n",
      "235/235 [==============================] - 28s 118ms/step - loss: 0.2011 - val_loss: 0.1998\n",
      "Epoch 33/50\n",
      "235/235 [==============================] - 30s 126ms/step - loss: 0.2010 - val_loss: 0.1998\n",
      "Epoch 34/50\n",
      "235/235 [==============================] - 28s 118ms/step - loss: 0.2009 - val_loss: 0.1997\n",
      "Epoch 35/50\n",
      "235/235 [==============================] - 25s 105ms/step - loss: 0.2008 - val_loss: 0.1996\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 25s 108ms/step - loss: 0.2007 - val_loss: 0.1996\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 24s 102ms/step - loss: 0.2006 - val_loss: 0.1994\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 25s 106ms/step - loss: 0.2006 - val_loss: 0.1994\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 24s 102ms/step - loss: 0.2005 - val_loss: 0.1993\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 24s 101ms/step - loss: 0.2004 - val_loss: 0.1993\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 20s 85ms/step - loss: 0.2003 - val_loss: 0.1992\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 18s 78ms/step - loss: 0.2003 - val_loss: 0.1991\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 19s 82ms/step - loss: 0.2002 - val_loss: 0.1991\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 19s 81ms/step - loss: 0.2001 - val_loss: 0.1989\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 22s 93ms/step - loss: 0.2001 - val_loss: 0.1989\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 20s 84ms/step - loss: 0.2000 - val_loss: 0.1988\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 20s 85ms/step - loss: 0.2000 - val_loss: 0.1988\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 20s 85ms/step - loss: 0.1999 - val_loss: 0.1986\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 20s 84ms/step - loss: 0.1998 - val_loss: 0.1987\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 18s 78ms/step - loss: 0.1998 - val_loss: 0.1986\n",
      "1595258110\n",
      "Wall time: 21min 15s\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+WH4yJAAANdUlEQVR4nO3bvYoU3RaA4aqjYiCIgWMggoKZ+IM63oTIB6Kh3oQXYqq5YCZi4i2YzCQmJqJ+iYEzgRiIKFLnAo7CNGdrv6PPEw7FYsNmVTcvPfOyLBMAAAAALf9Z9wEAAAAA+F+iDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAEDQwVUePnDgwHLo0KFfdRZ+4tu3b9P379/nEbOOHz++nDlzZsQoVrS9vb27LMvGiFnzPC8j5rC6ZVns4j43chfd43q8e/du2t3dHbKL3qdrNWwXjx07tpw8eXLEKFbw/v376ePHjz4X97mRn4tHjx5dTpw4MWIUK/jw4cP06dMnn4v73w93caVoc+jQocnL9Pd79+7dsFlnzpyZtra2hs1j7+Z5/nfdZ6DDLq7PyF10j+uxubm57iMwxrBdPHny5PTo0aNR49iju3fvDpvlfbo+Iz8XT5w4Md2/f3/UOPbo3r176z4CY/xwF/17FAAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEGiDQAAAECQaAMAAAAQJNoAAAAABIk2AAAAAEHzsix7fnhjY2O5devWLzwOP/LkyZNpZ2dnHjHryJEjy7lz50aMYkVbW1vby7Jsjph18eLF5fnz5yNGsYLr169PL1++HLKL8zzv/eXLaMN20T2uz7IsQ3bx3Llzy+PHj0eMYkWXL1+2i3+AUbs4z/Ny4MCBEaNY0ffv3+3iH2DULh4+fHg5derUiFGs6M2bNz/cRb+0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAg6uMrDp0+fnh4+fPirzsJPbG1tDZv1+fPnofNYj1evXk3Xrl1b9zH+Oru7u8NmnT9/fnr27Nmweezd2bNnh806f/789PTp02Hz2JubN28Om/X69evpxo0bw+axHlevXvX9Zg02NzeHzdrY2Jhu3749bB579+DBg2Gz7OJ6jNzFCxcuuMM1mef5h3/3SxsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIOjgKg9vb29P8zz/qrPwG1y5cmV68eLFuo/xVzp8+PCwWZcuXZq2traGzWNvNjc3h81almX6+vXrsHmsx9u3b6c7d+6s+xh/nbdv3w6bdfHiRe/TNRn5nfLly5fTqVOnhs1jbz58+DBs1s7OzvTgwYNh81iPL1++TK9evVr3Mf46X758GTbr9evX0z///DNsHv8/v7QBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAILmZVn2/vA870zT9O+vOw4/cXpZlo0Rg9zhWrnH/c8d/hnc4/7nDv8M7nH/c4d/Bve4/7nDP8MP73GlaAMAAADA7+HfowAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACC/gsAyDCuGGSCkwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dW7NVxdWGBzl8ORqjgkcQBEEUlUiUaKK5SuUilVTu8wfy73KwUolXJmo0lYjgAVAQEEQhihpNjDmH7yLl9OmXPdu1txuYa+3nueq1Z+916NGju+es8Y6x7vz58yUiIiIiIiIiItPiU5f7C4iIiIiIiIiIyIX40EZEREREREREZIL40EZEREREREREZIL40EZEREREREREZIL40EZEREREREREZIL40EZEREREREREZIJ8Zjmd161bZ33wy8T58+fXrcb7aMPLylvnz5/fsBpvpB0vH/riQqAvLgD64kKgLy4A+uJCoC8uAPriQrCkLxppI3LpOHW5v4CIVJW+KDIV9EWRaaAvikyDJX3RhzYiIiIiIiIiIhPEhzYiIiIiIiIiIhPEhzYiIiIiIiIiIhPEhzYiIiIiIiIiIhNkWdWjRJbDunXrRl+z/elPf3q036c+9amZ+uVn/fe//x3a//nPf5ZsV1X9+9//XvJ/ZHVJ+3zSfsn58+eXbMvFYVbfTuhj2ml1GVs3e9eW42+0F9fRXDe163RYjfVULg+z2k5bXXxm3e9Wup72zi/a99Iw6/7ZszFt1dsXtamsFCNtREREREREREQmiA9tREREREREREQmiPIo+UR85jPtFPr85z8/tK+++urm2oYNG4b2Nddcs+Tf8z0YipjSpn/+859D+/3332+uvffee0P7rbfeGtrvvPNO0+/Pf/7z0P773//eXOPnGc74EbOGBs8qyViN0NOerQxLXTmzShr5umfDnlRR6dTyGfOxDO8eu5b9SPoYX6/EZ2V59MZ11nW3d623LvZ8UZteHFYqZZPVY8x3ZpWbrnQ9de+7NMx6Ru2db3o27p1vlBTLamCkjYiIiIiIiIjIBPGhjYiIiIiIiIjIBPGhjYiIiIiIiIjIBDGnjcwEtZ9f+tKXhvbGjRubfvfcc8/Qvu+++5pr27dvH9qbNm0a2pn75v/+7/+GNrWkqQNlDhrmsKmqOnv27NB+9tlnh/aTTz7Z9Dt06NDQPnfuXHPtb3/729Du6VEXkV5ehDENd+Y3+tznPje0mafoC1/4wuj7JSzJzhxG//jHP5p+fP2vf/2ruTZmO3M1/I+edp++SBt+8YtfbPp99rOfHX1/jjl9iu2q1m6pBx/LZ7QWbLaSvFGpyadv5jVCf+vlYOiN+0pK4K4lX+yNSS8/Bn2MfpnrKfvlZ9GGvfWU13I9NRfHJ2PM/uaGujjMumbma66Tub/R//LcQ3o+xrWW9l2N3EbOl4/o2XvMxrRvVWvj3nziuaW3bvby+fEa50iVdp0SK80X90kx0kZEREREREREZIL40EZEREREREREZIIoj/oEzBq2TnqlNacMpS7XXnvt0L799tubfnv37h3au3fvbq7dfPPNQ3v9+vVDO8O7GbLYC+VnGOGVV145+h6vvfbaTJ81K2nbeQ1ZnDVMP/uN2SfDhPma8+fLX/5y04+vMwyZIcV/+ctfhjZLtVf1y4GP2afne/Nq0zF6NqRP3HTTTU2/nTt3Dm36+tatW5t+DCdO2/A1ffH06dNNP75+++23m2u0/axSuHliVjli+hj9hT6W8jXaZ0wiU9VKTvPamCym52+zyp5WKsWaF2Y9K/TWU9r3K1/5ytDmXlrVyo3zPWjfP/3pT0P73XffbfqlD5P0uQ+Z17PNajFrueix8sG9c8VyfGys36zXFqEccU8qyrWQPlXVrps8U371q19t+tH/+FkffPBB04/7GP2tquqvf/3r0OZam5KYla61Y/0WkVklcCl74tmHaR9o36p2XuQ9BOcQxznl3zzDvP/++801zgWu0WxXtfNkra2vF5Ox+dK7p8l7FV7jPMg1gefXldjQSBsRERERERERkQniQxsRERERERERkQky9/KoWSvbZIgkX/dkIisNTVxJJukphbtlSO9YhvwM32PY3zvvvNNcYzgZQwXTNhw7fm5+J4YKZnj3yZMnh/aRI0eG9ptvvtn0Y+jaWsjUPjbXe1UUetfG5nlVO55jIfVVrTyKIapV4yHEyUqqQi16tQ7aKUM5Kal44IEHhvb3vve9ph+v3XjjjUM75Tccf9qsqvV1houfOnWq6feHP/xhaO/bt6+5dvTo0aHNdSV9dp4Y87G0Fat2XXHFFc21G264YWjTPhnOz3WUa3TKYFiFL9f2sbDtDAPn6/TZsepUs1ZimBdWWrGmJ5cZk1FlKP/mzZuHdkpR33rrraFNKWFPbrrS8Z93G1bNboOq1k9pk6uuuqrpt2HDhiX7pVSHe2bKKSi16cmGeb6ZtapiMi9SRa5xlMHkmYLyQa6fVVW33HLL0N6yZcvQTvk9JTIcn5RC0N8oDa5q5cB//OMfh3aem/mePcnqrEzZhsth1qqX9LHrrruu6ce1kukb+PeqNiVErqmcd9wXWcG2qj3DHD9+vLn26quvDm3ey6Tce61VsZ2FsXmQ5yieWfMcxde0b0rhenB95TqcUn++zjPWLBhpIyIiIiIiIiIyQXxoIyIiIiIiIiIyQXxoIyIiIiIiIiIyQS55TptZywz3NIrUqKZe+Prrrx/aGzduHNrU+1e1GrbMi0BN6euvvz6033jjjaYfdaiZU2WsNFtqFKeqS8zvxTJl1Okm1N8ePny4uUZtKd8/dX1j+U9SX0gbpn6RGnDqRVNnSj14fu6YlnueNME9Tf5K8izk/81alpT2ueaaa5p+27ZtG9qZn4E5iDhPUidKv5oXH1ttemtmarmZq+aHP/zh0P7mN7/Z9ON6Snsy90lVu/5ljpOx8uK33nrr6PfNvDj02yydOi+kfegTY/kwqtoy7Cy7XlW1a9euoc09Lkubcj3kPpb7Fn0s/Yj2p2Y790Xmbjh37lxzjXbt5bmaF1aSI6yXh693PuI1zpf0beblyPwLnBdcQ/NzeSbqlSDmnOiVrJ46Y/td+hHPnjnu3Md27ty55N+r2vwYeW4h9I9cU3lG5fnm2LFjTT/m88uzD89q87hHZp5F2oY5vTJvzfbt24c218+q1lbMY5NzmfZgm/lyqtq8OHfffXdz7cyZM0P7hRdeWLJd1dp3Vl9cFMbOq1Wt72TeIp5bduzYMbTvueeept9dd901tLl/5h7Mz8p5x1xU3DPTZ3neefTRR5tr9MWxubXWoL25Duf9A32Y9uU+WNXOkZxLPKPys3L/pI/x/r+qXWtPnDgxtPM+mGcx3lfn+49hpI2IiIiIiIiIyATxoY2IiIiIiIiIyARZNXlUr7w24bUMM2NJrixZunXr1qHNcCiGt1W1oY8Mo8pwZYYlZQjaWOh3hp4eOHBgaB88eLC5xrBwhszNawgjwzI5XilToSwsw3EJxz/lUfwshiVmyCJLZ6ZMjmPO70g5VFUbgrycku5riV7oNMelJ/9gaCJLK1ZV3XbbbUObYf9V7RpBqVSGN/bkUbOW/O5JE+aBHH+upwzTrmrlUbt37x7aLAVe1foiZS/PPfdc04/hoDlfGKrMEOT0Wa7XOUcYptqT5E2ZtA/DqhmKn7/93nvvHdoPPfRQc41yKUrPUpbEvWtMrpavM+ScclSu31nKlv+XvteT3cw7vfWjt39wvZq1TDb9IeUflNOlvIfyNK61uWbSNrOup/NELzye45LnUJYCTqkF19Q77rhjaKdkhvA8kmcT+lh+D9qc7SxNzfUxy4bz/WnvtO9U7d0ruc5xSP+gBIr3C1Wt73DeU8pU1aZNYL+0NVM05DXOJe7VOQ94fk3Z8BjzYsOPg37ak8Nl+gveF95///1De+/evU0/zgW+f8q/KUfMMwf3RbbTZ7lX53zi/sz9dFHsOEavRDfv73ju+dGPftT0u++++4Y2zy95X8m0HrRnVWtT7q25nvLMRt+uateBU6dODe3VtqGRNiIiIiIiIiIiE8SHNiIiIiIiIiIiE8SHNiIiIiIiIiIiE2TVctqMlaOsGtcLZz/mNMgyiSxFS11i6sr4Pahhy9Km1PemRnGsvFhq8KlrTb3lWJnvedWH83vyN6RukNeYGyivsZ3jT50jtb6pQ6d+kVrD/F7Ms5M2XIR8Qz1m/U1pg1nz2PTyONAnmL+EWu6qNl9Vfi61p5wX+bt682mR4finxp85oFK7f+eddw5trqGpK6a++rHHHhvajz/+eNPv9OnTo+/Bz2bZ016Olyy1yGvznnvoQzhO9A+WC66qevDBB4c282ZUtbpvlkLPMpPUzDNPW67f9NnMwcAcAtR657zjepv7LvM1sN8i7Iu9v8+aq6ZX/pvzhfkSmDOqqs3hkfsd7Tu2p+frRc+rUNWOO88VWV6WZXszpw1fc4/jPK+qOn78+NB+6aWXhjbPq1WtDTKfH8uN84yU+Vvow3nenvd1NPd5jhftmWdDnu9zLvO8QXscOnSo6cd8jRzXTZs2Nf3G1syqNv8JfTj9mfbt5QtdRMbKPle1PsFcRFVtnkTmPOG5p6odzzG/rGrtnfOO+XR433rttdc2/fh/vXtkMq/7IqENe2e+nPc/+MEPhvaPf/zjoc0S7lXt2PGek/fnVW0exizXzXWeuYfyDNS7B+FZh5/dy/u6EnsaaSMiIiIiIiIiMkF8aCMiIiIiIiIiMkEuijwqwy4ZAsXw3AwJY/hvhqUyBI3vkWW4GcbIcLcM02ZYEkNNq9pyjQy7y5B9vkeWDR+T5CxCqDHD/Fgyu+rCcL4xerIOhicz/D8lHnydYXcnT54c2rRFhirP4/h/EmYN5++Ny5iv90oy9so5U56ToYRjUryeLG9Wm/YkDPNIjj/HnGW9q1oJDsOvU9J49OjRof38888P7QwfZih5lrrk2k275drBdaBXvjx/57zC38gx68kHU/7Aef/yyy8P7aeeeqrpR3uxnGn6Eb9Tfg+GCjM0Pe3NfnmNeyilUvPue1X9dWcl+0zKDDnm9O20U2895d79wQcfDO2LsZ5Omd5+19vTeJbN0HnObZ43uW5WVf3ud78b2iwNm2WGuR5yDahqz0iUYXCdrLpwDhH+5rH2lEmZCkuYc205d+5c049jnvsMx/zVV18d2q+88krTj6W3KXtKKRbnSMplKPfh56bN8juSRVg3e/RKQnOu8z6hqr034HqY53/add++fUP74MGDTT/KaXqSGd5LUvJc1abkyPtR3j/2ZMPzIv0fk0SlDblm8r67qur73//+0KYNc0zefPPNof3EE08M7V/84hdNvxdffHFopzyN6VcogU3/4hqT8qtnnnlmaPO5BGXrVReee5eLkTYiIiIiIiIiIhPEhzYiIiIiIiIiIhPkosSZZ8UChhSxnaFqDPnMEE2GIx44cGBoMySpqmr//v1DmyH7GZLEsKw9e/Y015hhnCHJGb7K909ZAcdgXkLaVkKvik+GljFMjmHHKTtjqOOuXbuG9re//e2mH2UXKZOjrINh4Dk35yUU+GKwGr+dNs2QQ8qjGCq6ZcuWph9D/Zmlv6qt5sAwwwznXw0fm8e50KuuQP+46667mmu0B30i5RQMEWfVoQzvpd+nJIYhw1kFhcy6Zs5rSHj+Js43+g6riiz1mtBfuBdmeDclAQzTTj/i98jwbq7ZtHGGK9Pvc06OSQLkf3A95ThWtb5D32allKrW37h+VrVnGF67GOvpPEFf5FmRZ4eqdpx61UNYTY9V96raMyrX27QB/T5lkWPrRYbi987bYxUX52UfzDnK38d5zgqISa5dlDdREpP7XdrqQ/Isy/uW/B/Os7H5V9W3zbxL3D6OXtoNSm1SHkUJPscv7xMoVaRfnjhxounHs0lvf+YanXOL94sprRlbBxbhfqVnQ0rceK9X1VaToh+lVPHhhx8e2j/72c+Gds+GKSmmxJH3I/l9aaff//73zbWnn356aFOSudo2NNJGRERERERERGSC+NBGRERERERERGSC+NBGRERERERERGSCrFpOm5WUgc6/U/uVOjDqUqkHzJw21BKzjFqWzaP2NEt+s7wYNa6Zb4N5AlL7PGvp63lnOWWTx/LY5Piz5Nr9998/tO+9996mH+fI4cOHm2vM29ArsTav+TEuJ2N5VJiTqqrV3VNzfMsttzT9+B6z5o1Km86j1nc14LqW5V6ZmyvLrPP/uHbl+I+V/s1cJdRy83Or2lKOLHua5VG5XtN/q9r9ZVHWVo4nf1OOCzXzOe/PnDkztLlHZplb5rqgPjz9hj6bpU2Zi2rTpk1DO23FPbiXt2Ye82hcDLie0vaZH4M5hnbs2DG0t23b1vTjmNPWVW15VM6JzHey1uBc5FjkuY62ypwqHNvjx48P7czBQN/slb3nep7rN3MyMPcRc9hUtXsm19eq+c+7mGsGfw/tlmd/jhfzV+S13piM3e9kSWPm/socVXydpeVnZRHz2IyRduR6uHXr1ubaWC7SI0eONP34mvtW7mmcJ5lDjOcd5hrLPEjMkcTPyr7MaTOPftkjbbh+/fqhnTakXzE/zW9+85um389//vOhzZxFmUOKe2veg3zta18b2jxv5XmYZeF/+ctfNtd4/vqkZb17GGkjIiIiIiIiIjJBfGgjIiIiIiIiIjJBVk0eNRbqna8ZApqh0wxvzLBeltpieGl+FsP2GaqYZdruueeeof31r3+9ucYQN4a8pgSHcqmUc8mFoXC0zViJtao2tJH9MnyYc4TliKvasGOGnqb0oBeiPybDWPQw1I+DdqVNs+Q3w7sZzs+QyKo2jJFyqKrW13vlVtcqtEWW2t6+ffvQztBsjmVPdkZ/4ftzjaxqS9KmPIqyGs6JDOWnP6ev51oyj+Sc5drD/SNlhvSxHLMxGWiG23Mv5PulzI2ypyzBSZkb1+wM9SZpN34v7v9rWabKMeIamj5LOcDtt98+tPNsQxkMzy9VrZyOEpK1vp7y9/dkQ3yd0inC80f6GM80XF/TjgzZv/vuu5trlEcxhD9lHWtJUjxWNjvXTL7OMaGtKE/MflyjeZ7ZvXt3049+mmkAaG+eX3O+cM3M9XSsnPI823bsN+XZnWcOprSoaseJ947pH4R+mWcp7n3f+MY3mmuUS3FPy32RNk7pFOfkIsuG04Ys653yftqNY5lyU67XfI/0I66Z3/3ud5trlEvxbPzcc881/R555JGhneXjL5XEeP5PwiIiIiIiIiIiC4gPbUREREREREREJshFkUfNSoaeMgSR4aVVbagwQxMzLI7VSfg/DI2qqtq7d+/QzqpEDIvcv3//0E55FL/jomX5Xg1S/ka7MYwtQ+PHKsXknGC27ldffbW5xhA3fm7KBvg6JXkMd+tlA1+0EMYk7UO7Mpw/ZR0MX925c+doP4ZwpzyK4aw9CeKihAYvF85fhvdWteHYOSYMqX/nnXeW/HtV639cT3Pd5WuGvFaN+1/KD3ktq3AwrLZXkWieGJNHLWcvYQgw7Z3VZtiPa2NWiKIk6r777muu0cY9+QHXzbQV91buARm2fjGrL0wNjlFPbsqzTUotCOVRlENVtWutku6P4PrYq4RK38k5y/eg7W666aamX8qgPiRtSqlNyqPot6xKw4pgVX0J3KJJa8ZsmPO8JyOlbSjTT7kM/48SqJTOUG6a/swx5/tlP77OfZFrx6LIasZSaKRclGPLtbFHvgcr1XLPpBS1qj2/pvybc4ZyVMp7qlpJVEppxnwxz97zaFdK1fLsT7/KMWFqBO5bec7lmYX3HCmPeuCBB4b2gw8+2FzjWeSll14a2o8//njTj3Kp3vp/MTHSRkRERERERERkgvjQRkRERERERERkgvjQRkRERERERERkgqxaThsyq7YrcxpQI0adblWrhaMuMfMnMPcBr1G7WFW1Z8+eoU1tZFWbV+PIkSND+/XXX2/6XaoSX/MENZipv6UNqTfs5b6h1pe6xqqqQ4cODe1eThvOF+qU81rqwTkfqQ1PjTTne879edSgfhxjJdRT980y39T157xgriLmV6lqx32sBPtaY8zHUuvLeZr6ao4r19qzZ882/VimkrlLUmvN15mThf9Hf+iVpc5rY+Wh56lUdH5XzmeuV5lXiLbrrV/cxzIXx4YNG4Y2116W7q5qS19mXhzalXrznFu0d+rK+X1p79wDFiHHxhi9cvZj+W2q2pwnzHOR6yLPL/TfqjZvm3n4lobjkmc8jmeeOW688cahTRvnukwfYI6nPMtu27Zt9Br9m+s3/bKqzW+UzNPaOQtj57Cc57RpnuXoV1xrs9Q2c54wJyZzn1Rd6MNkrFR7+jPXySyZzLWW/zfP59CxXD+596VfkbG1Ld+D+Wjoi9wvq9q1N/Pd0Ha8R+zlZ+ydZRfNL0n6Ec/+p06daq7RVlzHMucT113mNso186GHHhraec/Pe0vmrn3++eebfjzr9O4DLyZG2oiIiIiIiIiITBAf2oiIiIiIiIiITJCLIo/qwRCinjwqpUgMRWLoGkPaqtqQOb7/+vXrm34sDZbhVgwxPXjw4NDOMHB+p3kKP1xNMtyN4aAZvkgbMMQwbcOwNoaGZjjys88+O7QzFJFhpJwjKQcgOZdYvo+/M0OOGYqZ5W8XgZ6NGeqdJUu3b98+tDkXMnSVkpy045g8aq36W9W4PCpDpzl2uXbxNcNSU5ZKOQBtk6GhDGVNiduYJCbnFdfhfH9+X4a3z9M8yO/K38i978SJE00/hutm+WCGF5Ncy8ZkNxkuzv9LaciYjC4lGT3ZMGUdynP+B23Tkztyfe3tRzw7cQ/Lvo7/R4xJa3I/57xnadiq1oe5Fqff03Zcy/M8zDK0aSuuy6dPn17y+1W1criUTK4V++f4cxzSd/ia+1aWKqZvcj3lnpuvc12kb3Jvpc2qxuVCVeNzqSexmfqeye/OsaWUqaodixx3jjX3yPQx+inbOc58nfbhvtgr+U0JTpaL5rxbhDQAY/LE/N3Hjx8f2nnmo705/nlu5P7J8wtlU1WtJCrPyjyzUhLFtbVqGvunkTYiIiIiIiIiIhPEhzYiIiIiIiIiIhPEhzYiIiIiIiIiIhPkkue0Ib1SfKkHpJaM+k/qfqvassNjZcKr2vwJqbPbt2/f0KZuObWMi6wJ7mliqSlMbSBz1WzevLm5Rk3h1VdfPbQzFwpfU6vKHDZVbY6b1ENyHrAM3JYtW5p+nAdpX77/yZMnh3aWBh+bm1XT1w+PQftnOV7mXaDudOPGjU0/akqpTc7S7cypkn7PUsj0t4tR0nJeygyP+WbmLHj33XeHdv4ejjN1u6nhpV6but/MA8Y8VPlZXCO4Xue8YgnUzM3w2muvDe2xEuJTJ78rteu0x/79+0ffI0tVci2mdj/3NO6tvdxBzHGT+VAI50Kum/TZzDVAG/d0/PNk1+WSmnzq9ekf3Leq2j2TvpN5jbhv5X5E3+mtp2uZXo46+kSuw8ztRDv27M18fpnDiP6R++Irr7wytGnvLPHOvTZzqiyy/Xslv7k+5RrH9Y/7Kc+TCfe39DeelXKMuYaynWsm1+icc3zNeZa/eZ7sO5abJ/cInsOPHTvWXOO49/KQsB9tnD47Ni+q2r2WtuP5q6o9S62l3I38PXku4dk/x4TnlF6+Id7n33bbbUM7S36P5butanPr8J4/7/WmkG/ISBsRERERERERkQniQxsRERERERERkQlyWeVRCUPBeiFoDHdLeQ5DQBnilvIohkgeOXKkufbEE08M7TNnzox+p0Vj1tKCDFtj+fWqqm3btg3trVu3NtduvvnmJf8vw4IZDnro0KGhzRC2qtbWWeKW0pwdO3YM7Z07dzb9GBKZUgGGHfP7MjQ5+2W543mV0I2F3VaNlwzOcERK5TgOKX3h+KXUguGItM9qhCn2JIBTZkzGlaH8DM9NG3KeMqw/yzezH9fglEfRh3Me8DVlA1lulaH9LLuY34PzYJ7Dh8fChnOd4151ww03NNeyFO2H5NxmGDj3wrQj/S/tw/egP+e6SblOzifOSf6ufI95tutS9Mrxcq+lbSiHqmr3OJ5FUjpDiWOGo/P/OObzLKdYbfjbc17SP1KKxGs8l9Jvqlob0/Z5vqRkgufQfE3ZRfob7Z/vv1ZsnGeFnpydZxueL3Mt5FjSd/Is21trad/eWsj3750nZy35PU/0ykW//vrrQzvnMn2O45mSbMr7uZdmP56fenJHzos853JdTtnN2HxaBB+ddT3N8+tYme+8N+U+yfu5m266qelHn8h18sCBA0ObqTDyO03BHkbaiIiIiIiIiIhMEB/aiIiIiIiIiIhMkEnJo0iv0gavZRggw9quueaaJf+/qg1B/9WvftVce+aZZ4Y2Q/KmEBp1MenJoxjayXC07du3N/327Nkzeo0yCYaRZsgcw9MYxsYQ0qpWmpOVNphF/I477hj9TvweGT7MEGd+x+yXIbbzQi+Etlc9iiHEDEekHCr/jyH8Kf9gmGtWlmJYMsOVVyo764UQz6N/8zekHzGktGdDtnv9uJ6m9JE+lhXaKFUkGaL69NNPD+2DBw8217gOz6vkMBmTA7PCUlX7e1M+SP/jesiw76o2TJ/tlFdxPqUkYKw6SYatUzZC6UbVuDxqClUZLhU9H6PdUvLL+cL1lBU4qlrpTFaWok3HQvLlI3rn0PQPjiHPS+lj9D/aOyX89KuscjQmaU3ZxVi1sKr53O9mpVc9ivM+/YPjR7ulL45Jp+jLVeOS0qpxmW/245zrVdkbay/1ekr0zqH8vbn38ayY9uHrMflpVXtPcssttwztPLNwLvSkctzvWIUor+V7LIrk++Po/bY8v9L2uWcSnnt27949tNPWPHuwQnRV1W9/+9uhTSn+FPdFI21ERERERERERCaID21ERERERERERCaID21ERERERERERCbIZHPaJGOaz9RDjpVTZJ6UqrbM8K9//evmGnNspM5ukenltKGmkzps5lSoajWimffiuuuuG9rMj5Blssd02JkzhTk2Nm7c2Fxj6XGWfstcD9Qc97Sq1FemzpGlG6emR5219GP2m7VELTXc+duZd4Ga/BdeeB9AbksAAArZSURBVKHpR9+k7rSqtcmsut9Zc/UsApyL1NnntSwfTD+gPTP/AvMqUDu8a9eupt+99947tDdv3jz6Wcy/sH///qbfU089tWS/qsXUfI/taZkzi3kXMqcB/492THtzX2TZ9fQH+lvml+L3Ze6HLH1MX88cG1wrF82ms647vfWUuVByD+Y84Bzp5QjL/EhrNY/QatErXzt2fsqS3zzH3HDDDUM7117m8Eg7cm+l/9G/qqaZk+Fy08tLRB/jWKav0Na9ctC98s1ja0LPZr1cNb2cNvPEmI+lrbg/bdiwobnG/Y75LrMMNPPv8d4lcxPxs/Nswtw1Tz755NA+duxY0683nxbFdstlJb8118mdO3cObd4H5vmF++LDDz/cXDt06NDQznP01DDSRkRERERERERkgvjQRkRERERERERkgsyNPIrMKuNhaPZzzz3X9Hv++eeHdoYXMxRuLYWqkZ5chmOSZbgZkpYhgHzN/8tyliyhzfDIlGLR1lnym+W6+Z0y9I3vn9IcllI9ffr0ku2qNlR5nkLOx0J883XOBf5Gji0lh1WttIbh3QxFrKp65ZVXluxXNS6VWw2Z07z69lj4cJZepuwwQ0oZRsoQ/R07djT9OA8oj+L/VLXSxyy/ye/Fst4/+clPmn4vv/zy0O5JvRaRWUvUpvxhzIdTksH35/6WZTHpz5SfVrXrI9fsLOtNH6b/VrXztScdmBfGJA65PtE2eWZhOVP2S1tTBsM5Qb+pqnrttdeGdq6n9KtF96mLQU/GMLY/UfJWVbV+/fqhzTNN+iz9NMsd0+fol/N0/pgCOV4cZ45rSvi5x3FvzX2Wtk/70hd7aRh657RFl0f15N9c29I/ONaUPWUaBV7jOswzfVXViRMnhvbhw4eba5REHThwYMnvV9Wu2fNUkv1ywXlPudqmTZuafpTq83/yvv6RRx4Z2o899lhzLW01ZYy0ERERERERERGZID60ERERERERERGZID60ERERERERERGZIHOT02Ysj01q8lnOlLkUsvTomTNnhnZqx6kvHMvlsoj0ys4yBw3HkvlI8v+Ym6bqwjwYH5IlY1majbriXp6AzKvA18xvk7+LORcyJwh/J0tB8vtVzZcecizvQmqlx8p6J9T+Hjx4sLnG9+RY0vfyPTIHxmrnXej58Lz491hOmyxx+OKLLw5t5pypqtq9e/fQZj6oXm6onj6fdspSl9R8//SnPx3a+/bta/r1SlsvIivJy5T+QH/hWpk5EsbyGzFPUVWbx4b2qGptzrk2a1nvqvnxsZVAezI/QlU7dnlmyRxQH5I51rh3cc9Jf2O+ocw5N+t62pubi2zDj2PWPZNrZdqXr7m39nIAZo6NsfxGve+7luj97l7uqbHcX3mWpa35Hrz/qLowDyPh+sqzZ5a2HssDlq8XJacN4e/LvYT7TK6V3J9694FcH7lHModNVVvWO3MyMpcjv0fvLLso9lkuvVypuWfSlzZv3jy09+7d2/RjGXf66f79+5t+jz766JL9quYrv5uRNiIiIiIiIiIiE8SHNiIiIiIiIiIiE2RS8qhe6CnD9K+66qqhvX379qbfzp07hzZDEzN8juGrWZKR4YkMm8oQqkULcePvyfJ6vMawvwyhP3v27NBOuQxt2JOdMdyXY54yHc6RDJ0ck3plKDFf53swPJnXFqUccS9MmKGKGbZIezFMP8N6OU4MUc3yjGNlvZPV9rd59V9+7zHJRFXV0aNHh3aG3lPit3Xr1iXbVa2Ugz6QNmQ48bPPPttce+GFF5bsl2vHWpBEjdHzxR70F/pR7p+E49wrP53rIfvSdlkOl2Hrud7Oq88tl+XITXmN4fqUOVW16ynD/LMf19pZfWqt2GW5pB17pdt5jqQEMSUyV1xxxdCmv6Wv8MyakmLKP3r751qy69i6mTbs3WcQ+lvKDDnmvVLO3CdTvnzu3Lmhzf2Yf8//68kd59XWs0rWcy0bS9lQ1aZtoB15f5LXaINTp041/SiBSqkibcx5kd93EWz1SUl/45p59dVXN9duvfXWob1nz56hTWl/VXtfSbunjI2SqHk+axppIyIiIiIiIiIyQXxoIyIiIiIiIiIyQS6rPKpXYSGzsF9//fVDe9euXUP77rvvbvpt3LhxaDME6uTJk00/Vs3Iz2J4McPnMqRtkUPceqF9HJPVqE4w63v0Mo/3bNELS1zEjPtk1tDTDKumNCKrqzEcmKG7+Vl8j7F2VTvXeuHdi2iflTBWPSqrnzFUN0OuDx8+PLQZvr9+/fqmH9dGhiPn+zFkOKWo/D9+X+35ET15VG+do+/0qsiMyUVT6kl/7sl4uCakLG9s/8zvO6+yUjK2PuUaRxumTIL2SN8Z60efSskEx7y3j8vSzOqLPdkNySoylBPSHmkrrqkp6+Br+l9Pur3o6y1/X0/2RNIf6GOUe6f8hrai/2VVGsqLe9IpVjvNvbUnQe9VlppXxlJS5NzmPpMVXSnfpbSJUpqq8eqLuaeNSaCq2rnA9qL726zQFzMVCe/rb7/99ubat771raF95513Dm1WuKxqbT12b1J14Z48rxhpIyIiIiIiIiIyQXxoIyIiIiIiIiIyQXxoIyIiIiIiIiIyQS5rTpvU/FETTC1oVVva+zvf+c7Qzpw2V1555dCmvjRL1DJXA8szVrXaferxUnM8a06VRcDcIovBWAnFtCm1xKkPp26U13r5gmbNW9P7vnIhY+W/83XmJaJGn+UtZ80vpV0+Ob1ypqRnE16jvVN3z7wL9NnMK9crTc3P4vsvJ6fKIs+bsVwM+TpzMzBPRZaRHnuP3nq6lvKYXAx6fsnXmSOBPsH1NvMUnT59emjTx3Je0K96c2YsZ1jV2rX/SvO7jK1xmUuGeTR4n5E5U5inM32b6zXzs+ReTduvhZLuY2eaPN/QJnlGfeONN5a8lnvpmK8vZ99aRBusJr3cmfSJG2+8sbm2Y8eOob1ly5ahnbbm+nrs2LGhnXls6VfznP/JSBsRERERERERkQniQxsRERERERERkQlyWeVRGaLEcLcsnUcYRnXdddc11xhuxRC5LP9FeUCvVG5PQmJYnCwK6Ys96d+soYWzSur0o+mgLabHSnwnw7spm5i1pHHvs5TgLI/emjlWjn1WHP/LQ447fY5n2ZTWsOT3rP6WzHN4/5QZO7OkPG1MepzSjV7pcUtFfzJ6clSZDpzPuRYePXp0aHNdrKo6e/bs0OY9/5tvvtn0O3z48NCm1D8l4oviV0baiIiIiIiIiIhMEB/aiIiIiIiIiIhMEB/aiIiIiIiIiIhMkHXL0XmtW7fukonCUuvLXDVf+cpXhvZVV13V9GPJvffee29oZxk9lgDPMolT1OufP39+tlq8H8OltKFcwDPnz5+/dzXeSDtePvTFhUBfXAD0xYVAX1wA9MWFQF9cAPTFhWBJXzTSRkRERERERERkgvjQRkRERERERERkgiy35PdbVXXqY3utAilLYom9t99+e8n2ArN5Fd/rktlQLkA7zj/acDHQjvOPNlwMtOP8ow0XA+04/2jDxWBJOy4rp42IiIiIiIiIiFwalEeJiIiIiIiIiEwQH9qIiIiIiIiIiEwQH9qIiIiIiIiIiEwQH9qIiIiIiIiIiEwQH9qIiIiIiIiIiEwQH9qIiIiIiIiIiEwQH9qIiIiIiIiIiEwQH9qIiIiIiIiIiEwQH9qIiIiIiIiIiEyQ/wdzHvZxOBBFugAAAABJRU5ErkJggg==\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: 11.2 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": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0339, PSNR (dB) = 14.70\n",
      "MSE = 0.0578, PSNR (dB) = 12.38\n",
      "MSE = 0.0160, PSNR (dB) = 17.95\n",
      "MSE = 0.0490, PSNR (dB) = 13.10\n",
      "MSE = 0.0386, PSNR (dB) = 14.14\n",
      "MSE = 0.0156, PSNR (dB) = 18.08\n",
      "MSE = 0.0516, PSNR (dB) = 12.88\n",
      "MSE = 0.0546, PSNR (dB) = 12.63\n",
      "MSE = 0.0813, PSNR (dB) = 10.90\n",
      "MSE = 0.0352, PSNR (dB) = 14.54\n",
      "Wall time: 2.99 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": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4573\n",
      "0.2892\n",
      "0.4081\n",
      "0.4795\n",
      "0.2979\n",
      "0.5250\n",
      "0.2104\n",
      "0.1610\n",
      "0.1476\n",
      "0.5952\n",
      "Wall time: 202 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": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7183\n",
      "0.5042\n",
      "0.9115\n",
      "0.6390\n",
      "0.5384\n",
      "0.9322\n",
      "0.5599\n",
      "0.3219\n",
      "0.1908\n",
      "0.8264\n",
      "Wall time: 155 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/tensorflow/models/blob/master/research/compression/image_encoder/msssim.py\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from scipy.ndimage.filters import convolve\n",
    "\n",
    "def _FSpecialGauss(size, sigma):\n",
    "  \"\"\"Function to mimic the 'fspecial' gaussian MATLAB function.\"\"\"\n",
    "  radius = size // 2\n",
    "  offset = 0.0\n",
    "  start, stop = -radius, radius + 1\n",
    "  if size % 2 == 0:\n",
    "    offset = 0.5\n",
    "    stop -= 1\n",
    "  x, y = np.mgrid[offset + start:stop, offset + start:stop]\n",
    "  assert len(x) == size\n",
    "  g = np.exp(-((x**2 + y**2)/(2.0 * sigma**2)))\n",
    "  return g / g.sum()\n",
    "\n",
    "\n",
    "def _SSIMForMultiScale(img1, img2, max_val=255, filter_size=11,\n",
    "                       filter_sigma=1.5, k1=0.01, k2=0.03):\n",
    "  \"\"\"Return the Structural Similarity Map between `img1` and `img2`.\n",
    "  This function attempts to match the functionality of ssim_index_new.m by\n",
    "  Zhou Wang: http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "  Returns:\n",
    "    Pair containing the mean SSIM and contrast sensitivity between `img1` and\n",
    "    `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  img1 = img1.astype(np.float64)\n",
    "  img2 = img2.astype(np.float64)\n",
    "  _, height, width, _ = img1.shape\n",
    "\n",
    "  # Filter size can't be larger than height or width of images.\n",
    "  size = min(filter_size, height, width)\n",
    "\n",
    "  # Scale down sigma if a smaller filter size is used.\n",
    "  sigma = size * filter_sigma / filter_size if filter_size else 0\n",
    "\n",
    "  if filter_size:\n",
    "    window = np.reshape(_FSpecialGauss(size, sigma), (1, size, size, 1))\n",
    "    mu1 = signal.fftconvolve(img1, window, mode='valid')\n",
    "    mu2 = signal.fftconvolve(img2, window, mode='valid')\n",
    "    sigma11 = signal.fftconvolve(img1 * img1, window, mode='valid')\n",
    "    sigma22 = signal.fftconvolve(img2 * img2, window, mode='valid')\n",
    "    sigma12 = signal.fftconvolve(img1 * img2, window, mode='valid')\n",
    "  else:\n",
    "    # Empty blur kernel so no need to convolve.\n",
    "    mu1, mu2 = img1, img2\n",
    "    sigma11 = img1 * img1\n",
    "    sigma22 = img2 * img2\n",
    "    sigma12 = img1 * img2\n",
    "\n",
    "  mu11 = mu1 * mu1\n",
    "  mu22 = mu2 * mu2\n",
    "  mu12 = mu1 * mu2\n",
    "  sigma11 -= mu11\n",
    "  sigma22 -= mu22\n",
    "  sigma12 -= mu12\n",
    "\n",
    "  # Calculate intermediate values used by both ssim and cs_map.\n",
    "  c1 = (k1 * max_val) ** 2\n",
    "  c2 = (k2 * max_val) ** 2\n",
    "  v1 = 2.0 * sigma12 + c2\n",
    "  v2 = sigma11 + sigma22 + c2\n",
    "  ssim = np.mean((((2.0 * mu12 + c1) * v1) / ((mu11 + mu22 + c1) * v2)))\n",
    "  cs = np.mean(v1 / v2)\n",
    "  return ssim, cs\n",
    "\n",
    "\n",
    "def MultiScaleSSIM(img1, img2, max_val=255, filter_size=11, filter_sigma=1.5,\n",
    "                   k1=0.01, k2=0.03, weights=None):\n",
    "  \"\"\"Return the MS-SSIM score between `img1` and `img2`.\n",
    "  This function implements Multi-Scale Structural Similarity (MS-SSIM) Image\n",
    "  Quality Assessment according to Zhou Wang's paper, \"Multi-scale structural\n",
    "  similarity for image quality assessment\" (2003).\n",
    "  Link: https://ece.uwaterloo.ca/~z70wang/publications/msssim.pdf\n",
    "  Author's MATLAB implementation:\n",
    "  http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "    weights: List of weights for each level; if none, use five levels and the\n",
    "      weights from the original paper.\n",
    "  Returns:\n",
    "    MS-SSIM score between `img1` and `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  # Note: default weights don't sum to 1.0 but do match the paper / matlab code.\n",
    "  weights = np.array(weights if weights else\n",
    "                     [0.0448, 0.2856, 0.3001, 0.2363, 0.1333])\n",
    "  levels = weights.size\n",
    "  downsample_filter = np.ones((1, 2, 2, 1)) / 4.0\n",
    "  im1, im2 = [x.astype(np.float64) for x in [img1, img2]]\n",
    "  mssim = np.array([])\n",
    "  mcs = np.array([])\n",
    "  for _ in range(levels):\n",
    "    ssim, cs = _SSIMForMultiScale(\n",
    "        im1, im2, max_val=max_val, filter_size=filter_size,\n",
    "        filter_sigma=filter_sigma, k1=k1, k2=k2)\n",
    "    mssim = np.append(mssim, ssim)\n",
    "    mcs = np.append(mcs, cs)\n",
    "    filtered = [convolve(im, downsample_filter, mode='reflect')\n",
    "                for im in [im1, im2]]\n",
    "    im1, im2 = [x[:, ::2, ::2, :] for x in filtered]\n",
    "  return (np.prod(mcs[0:levels-1] ** weights[0:levels-1]) *\n",
    "          (mssim[levels-1] ** weights[levels-1]))\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = MultiScaleSSIM(original, contrast, max_val=1.0)\n",
    "    print(f'{a:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Simple Autoencoders in Keras - 784-4-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
}