Simple Autoencoders in Keras - 784-8-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 = 8"
   ]
  },
  {
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1595257560\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 1:47 - loss: 0.6925WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.417997). Check your callbacks.\n",
      "235/235 [==============================] - 6s 25ms/step - loss: 0.3335 - val_loss: 0.2523\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 5s 23ms/step - loss: 0.2275 - val_loss: 0.2056\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 6s 25ms/step - loss: 0.1978 - val_loss: 0.1895\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 5s 21ms/step - loss: 0.1870 - val_loss: 0.1826\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 5s 21ms/step - loss: 0.1822 - val_loss: 0.1793\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 5s 22ms/step - loss: 0.1796 - val_loss: 0.1773\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 5s 23ms/step - loss: 0.1777 - val_loss: 0.1757\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 6s 25ms/step - loss: 0.1762 - val_loss: 0.1743\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 6s 27ms/step - loss: 0.1749 - val_loss: 0.1732\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 6s 25ms/step - loss: 0.1739 - val_loss: 0.1722\n",
      "Epoch 11/50\n",
      "235/235 [==============================] - 6s 25ms/step - loss: 0.1729 - val_loss: 0.1712\n",
      "Epoch 12/50\n",
      "235/235 [==============================] - 6s 24ms/step - loss: 0.1720 - val_loss: 0.1703\n",
      "Epoch 13/50\n",
      "235/235 [==============================] - 6s 24ms/step - loss: 0.1713 - val_loss: 0.1696\n",
      "Epoch 14/50\n",
      "235/235 [==============================] - 8s 33ms/step - loss: 0.1706 - val_loss: 0.1690\n",
      "Epoch 15/50\n",
      "235/235 [==============================] - 8s 35ms/step - loss: 0.1700 - val_loss: 0.1684\n",
      "Epoch 16/50\n",
      "235/235 [==============================] - 11s 48ms/step - loss: 0.1695 - val_loss: 0.1679\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 14s 62ms/step - loss: 0.1690 - val_loss: 0.1674\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 14s 58ms/step - loss: 0.1686 - val_loss: 0.1670\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 17s 74ms/step - loss: 0.1683 - val_loss: 0.1667\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 14s 59ms/step - loss: 0.1680 - val_loss: 0.1664\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 16s 67ms/step - loss: 0.1677 - val_loss: 0.1662\n",
      "Epoch 22/50\n",
      "235/235 [==============================] - 16s 69ms/step - loss: 0.1675 - val_loss: 0.1660\n",
      "Epoch 23/50\n",
      "235/235 [==============================] - 17s 71ms/step - loss: 0.1673 - val_loss: 0.1657\n",
      "Epoch 24/50\n",
      "235/235 [==============================] - 14s 61ms/step - loss: 0.1671 - val_loss: 0.1655\n",
      "Epoch 25/50\n",
      "235/235 [==============================] - 17s 71ms/step - loss: 0.1669 - val_loss: 0.1654\n",
      "Epoch 26/50\n",
      "235/235 [==============================] - 20s 83ms/step - loss: 0.1668 - val_loss: 0.1654\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 16s 70ms/step - loss: 0.1666 - val_loss: 0.1651\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 16s 67ms/step - loss: 0.1665 - val_loss: 0.1650\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 16s 68ms/step - loss: 0.1664 - val_loss: 0.1649\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 14s 60ms/step - loss: 0.1662 - val_loss: 0.1648\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 18s 76ms/step - loss: 0.1661 - val_loss: 0.1647\n",
      "Epoch 32/50\n",
      "235/235 [==============================] - 24s 102ms/step - loss: 0.1660 - val_loss: 0.1645\n",
      "Epoch 33/50\n",
      "235/235 [==============================] - 26s 111ms/step - loss: 0.1659 - val_loss: 0.1644\n",
      "Epoch 34/50\n",
      "235/235 [==============================] - 17s 73ms/step - loss: 0.1658 - val_loss: 0.1643\n",
      "Epoch 35/50\n",
      "235/235 [==============================] - 17s 70ms/step - loss: 0.1657 - val_loss: 0.1643\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 17s 70ms/step - loss: 0.1656 - val_loss: 0.1642\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 19s 81ms/step - loss: 0.1655 - val_loss: 0.1639\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 21s 90ms/step - loss: 0.1654 - val_loss: 0.1639\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 19s 80ms/step - loss: 0.1653 - val_loss: 0.1639\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 18s 78ms/step - loss: 0.1653 - val_loss: 0.1638\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 19s 81ms/step - loss: 0.1652 - val_loss: 0.1638\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 26s 112ms/step - loss: 0.1651 - val_loss: 0.1638\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 18s 76ms/step - loss: 0.1651 - val_loss: 0.1636\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 21s 89ms/step - loss: 0.1650 - val_loss: 0.1636\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 23s 97ms/step - loss: 0.1650 - val_loss: 0.1635\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 22s 95ms/step - loss: 0.1649 - val_loss: 0.1635\n",
      "Epoch 47/50\n",
      "129/235 [===============>..............] - ETA: 7s - loss: 0.1648"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import os\n",
    "import time\n",
    "\n",
    "CURRENT_TIME = int(time.time())\n",
    "print(CURRENT_TIME)\n",
    "\n",
    "try:\n",
    "    os.makedirs(f'autoencoder_mnist/logs-{CURRENT_TIME}')\n",
    "except Exception:\n",
    "    pass\n",
    "\n",
    "# Run this in a terminal\n",
    "#tensorboard --logdir=autoencoder_mnist/logs\n",
    "\n",
    "CALLBACKS=[\n",
    "    TensorBoard(log_dir=f'autoencoder_mnist/logs-{CURRENT_TIME}')\n",
    "]\n",
    "\n",
    "# Fit data\n",
    "autoencoder.fit(x_train, x_train, epochs=50, batch_size=256, shuffle=True, validation_data=(x_test, x_test), callbacks=CALLBACKS)\n",
    "\n",
    "print(CURRENT_TIME)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoded_imgs = encoder.predict(x_test)\n",
    "decoded_imgs = decoder.predict(encoded_imgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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": "iVBORw0KGgoAAAANSUhEUgAABEUAAAHECAYAAADfz5RiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAARpElEQVR4nO3bsWud5QLH8Sf35FTMIfVA0qjB2IproSBZsoguIoguDg4OdskgSlUsuAiCq4KbVuggKkgnwTo5WFDpUFoFLRYxosFiESU1CbE1aXLuP3B7ueG8531u8/t81vcN/A5P3gS+vGdsMBgUAAAAgDT/qj0AAAAAoAZRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBrfzc2dTmcwPr6rH7kl3Lhxo2xvb4/V3tGGfr8/mJ2drT1jJC5duvTnYDA4UHtHG3q93qDf79ee0bi//vqrbGxsRDyL09PTg0OHDtWeMRIXLlyIeRbHxsYGY2N771d2MBiUwWCw9z7Yf3DbbbcNer1e7RkjcfXq1ahnsfaGUUl5FvfyGZZSYp7F22+/fTA5OVl7RuPW19fLtWvXIp7Fffv2DSYmJmrPaNzff/9dNjc3/+MZ7qpwjI+Pl3vuuaeZVf9HLl++XHtCa2ZnZ8sHH3xQe8ZIzM/PL9fe0JZ+v1+effbZ2jMa984779Se0JpDhw6V8+fP154xEmNjYzHP4tjYWOl2u7VnNG5ra6v2hNb0er3yyCOP1J4xEqdOnYp5FuH/XMyzODk5WZ566qnaMxp36tSp2hNaMzExUR588MHaMxr3xRdf3PSar88AAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBI47u5udvtlgMHDoxqSzW///577Qmt+fXXX8vx48drz2BInU6n9Pv92jMa1+l0ak9ozTfffFPuuOOO2jMYUqfT2ZPnePXq1doTWrNX/56m6fV65fDhw7VnNO7ixYu1J7TmwIED5cknn6w9YyROnDhRe0JrVldXy6efflp7RuNWV1drT2jN6upqOX36dO0ZrfKmCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACDS+G5u7vV6ZWFhYVRbqvn5559rT2jN5uZm1Ofdq2ZmZsrzzz9fe0bj3nvvvdoTWrOzs1PW1tZqz2BIN27cKH/88UftGQzh4MGD5cSJE7VnjMS7775be0Jrut1uufvuu2vPaNwPP/xQe0JrNjY2ytdff117BkPqdrtldna29ozGJf2vn5ubK8ePH689o3FvvvnmTa95UwQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBI47u5eW5urrz11luj2lLNl19+WXtCazY3N8vy8nLtGQxpaWmpPPHEE7VnNG5paan2hNYcPHiwvPrqq7VnjMTi4mLtCa2Zmpoqjz/+eO0ZjTt9+nTtCa25fPlyeeWVV2rPYEj3339/+fjjj2vPaNz8/HztCa2ZmpoqTz/9dO0ZI3Hu3LnaE1qzsbFRzp49W3sGQ5iZmSnHjh2rPaNx77///k2veVMEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASOO7uXl5ebksLi6Oaks1y8vLtSe0Zm5urrz88su1Z4zEiy++WHtCa7a3t8va2lrtGY3b3t6uPaE1W1tb5bfffqs9gyH1er2ysLBQe0bjPv/889oTWrOxsVHOnj1bewZDWllZKR9++GHtGY1bWVmpPaE1MzMz5dixY7VnjMQLL7xQe0Jr+v1+efjhh2vPaNyZM2dqT2jNlStXyuuvv157RuOuXLly02veFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEGt/NzX/++Wc5efLkqLbQgk6nU/r9fu0ZDGliYqIcOXKk9ozGXbp0qfaE1nS73TI3N1d7BkNaX18vZ86cqT2jcevr67UntKbX65WFhYXaM0biq6++qj2hNTs7O+X69eu1ZzRuZ2en9oTWXLt2rXz77be1ZzCkO++8s7z00ku1ZzTu+++/rz2hNdevXy8//vhj7RmN+2//I7wpAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQKTx3dw8MTFRDh8+PKot1Vy8eLH2hNZMTU2VZ555pvaMkTh69GjtCa2ZmJgoDzzwQO0Zjfvkk09qT2jNvn37yuzsbO0ZDGl6enpP/k29cOFC7QmtmZycLA899FDtGSPxxhtv1J7QmsFgULa2tmrPaNxgMKg9oTWdTqf0+/3aMxjSyspK+eijj2rPaNzKykrtCa1ZW1srn332We0ZjVtbW7vpNW+KAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKN7+bmmZmZ8txzz41qSzWvvfZa7Qmt+e6778p9991XewZD+uWXX8rRo0drz2AI3W633HXXXbVnMKT9+/eXRx99tPaMxu3fv7/2hNYsLS2Vxx57rPYMhjQ9PV0WFxdrz2jcyZMna09ozT///FN++umn2jMY0r333lvefvvt2jMad+7cudoTWnPkyJFy/vz52jMaNz8/f9Nr3hQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAgkigCAAAARBJFAAAAgEiiCAAAABBJFAEAAAAiiSIAAABAJFEEAAAAiCSKAAAAAJFEEQAAACCSKAIAAABEEkUAAACASKIIAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQSRQBAAAAIokiAAAAQCRRBAAAAIgkigAAAACRRBEAAAAg0thgMPjfbx4b+6OUsjy6OdUcHAwGB2qPaMMePsNSnONe4Az3Bud463OGe4NzvPU5w73BOd76nOGt76ZnuKsoAgAAALBX+PoMAAAAEEkUAQAAACKJIgAAAEAkUQQAAACIJIoAAAAAkUQRAAAAIJIoAgAAAEQSRQAAAIBIoggAAAAQ6d+FL7vBSdOWRwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2d2bMdVfXHV5xAQAYBE8QAScgACXMCFMggFIiUVVKWA+WT+qDlm3+Fb/4DVukLT5RVSiFQMiklASohISSBMCSEJAwJIQgi4ICa38Ov2Hz2yu2dcy/3Jn36fD5P69A75/bZq9fu3c36rjXv4MGDISIiIiIiIiIi/eJTR/sERERERERERETkUHxpIyIiIiIiIiLSQ3xpIyIiIiIiIiLSQ3xpIyIiIiIiIiLSQ3xpIyIiIiIiIiLSQ3xpIyIiIiIiIiLSQz4zncHz5s2zP/hR4uDBg/Nm43v04VHlwMGDB0+fjS/Sj0cPY3EQGIsDwFgcBMbiADAWB4GxOACMxUEwZSyaaSNy5Nh9tE9ARCLCWBTpC8aiSD8wFkX6wZSx6EsbEREREREREZEe4ksbEREREREREZEe4ksbEREREREREZEe4ksbEREREREREZEeMq3uUSLTYd68eZ2fu+yIiE9/+tPFPniwu3j5Zz7z8eX7v//9rzr23//+d8rvyN/Hf9f6WyIi40JeUz/1qY///0xrPXQNFBGpyetpF66f40vLx61jrfupyGxjpo2IiIiIiIiISA/xpY2IiIiIiIiISA9RHiXThvKlz33uc9Wx4447rtinn356dWzJkiVT2kuXLq3GnXTSScX+/Oc/X+wvfOEL1bhjjjmm2P/617+qY3//+9+L/dJLLxV7w4YN1bh169YVe//+/dWxf//738U27fFjulJFZ5peOuoxZW4i3bQkp4RSqTyOcZQlp8aViEwKM9mz5P+udGZ8oO/yc80JJ5xQbD53RET85z//Kfb7779f7PxMwnFeCzJTzLQREREREREREekhvrQREREREREREekhvrQREREREREREekh1rSRkWAdBNaWOeuss6pxF198cbGvv/76zmMLFiwoNmvYRER89rOfLTbberc0xh9++GH1+d13353yb7355pvVuKeffnqk759kWu2DeYy+iqh1wccee2yxWacof85aYrZuZ50i+jci4oMPPig2tcMRh9bmkJqWD+mb0047rdhnnHFGNY4xm2tD7d27t9jvvfdesbOf1HkfnhyL/My4pJ0/c54ZX/nYTM+phT4+lK5aRHleWUuO/uR/j6hjOB/j/POeyfptEXVsun7KkGjFWNcamuOIMD7yepo/z4RR203L4Wn5m3vU/ExyyimnFDvXtOHaydo3rG8TEfGPf/xjyn8TUa/FrbpyImbaiIiIiIiIiIj0EF/aiIiIiIiIiIj0EOVRMhJMI6SEZdTWshER//znP4vNdMDcGm8maalMPYyIeO2114q9du3aYj/yyCPVuH379hU7pyxOUuppS3bRSimlfyiRiajTTZk2+sUvfrEa9+Uvf7nYxx9/fHWMkqg9e/YUO/u71ca4q234JNPV3pJpwBERa9asKfY3vvGNYq9atarz+5555pnq2H333VfsDRs2FPutt96qxnFNmGQ/jRpvEbXvKGWjhDUi4rjjjis21+F33nmnGkeZYSu1vxVvLfkVP3P9Hrq/ZyLJyFJFrq/0e14zv/SlLxU7p/nzXvv6668X+8CBA9U4pvZn3wzdV6My11LFljRiVEkdyX7rir8h+LfLFxH1Gprl2NyzdNkRdWwyVrJsm2vtqFKX1vy3YnFoPpwteC3kPeqpp55a7BUrVhT7/PPPr8adeOKJxaZPI2qf837617/+tRrHz/m+S9k497a57INtw+ee2SiTMZe+MdNGRERERERERKSH+NJGRERERERERKSH9FYe1Uo9HeW/52M5rZxpcvkY01SZntbqSjP0VDX+VqbysTNMRD3nOR10+/btxT799NOLzXTufIyp39lP9EfuWLN+/fpi//nPf57yHCJqOcCkVWpvpRC3Yqfr3+UYoH8YUzklmSn8ubMUU0Vb19ZspAYPWUaV10mmd5988snFvummm6pxP/vZz4q9fPnyYud0ca6T7NYWUXdboN/WrVtXjWPKcF5rh05LTsF7VZY9UVq4cuXKYuf0bkoSKTPcuHFjNY7Hcho44fWTY5H/jutrRB3PXevDuDKqJKO1no66JvOaYFe3iIgLLrig2DkWmaLP72t145tkWtKmfB+jBJFravYPP/M+QylwRHtvws+Mo+w3SneypJhycMbfdGRafaIrdrLMkPejvJ7SbzyWJYj0Pef8b3/7WzWOx7L8nlJF2tlPrU5uQ9unzBZd8u9zzjmnGnf77bcX+8Ybbyx2vmZefvnlYr/44ovVMd7vKC9vPZvmkhBd+1z9e3i64j5L4ShxyyUaGN/0PZ91I+r7ZD7G+GbMZonbJ11PzbQREREREREREekhvrQREREREREREekhvrQREREREREREekhR6SmDTXc1ItlTTA/U3fa+ne5tgJ1a9T7L126tBp31llnFZta5Ihaq7Z79+5ib9mypRq3bdu2YueWmdS3UZc4rppUnjc1nFmnS80f5y6i9hX9m31NLTHr22Q9KnWJuf4CfcPzyHrhcZn/2aKrdsZ02mR31Y/J13ZXO/D58+dX4xiL+Tyo82/VrRkX3f3RIs8rY/Haa68t9s9//vNq3LJly4pNH7auify3GMMXXnhhsbP+/9lnny12rrExhJonma5aJrnmCe9Pixcvro6xBtHXvva1YnPOIyL27dtXbLZ6zn8rr8Vdx7hGZ+04fZdrjTFO6dNW2/A+01U7I9dEaLVl5px02fnfteo0rFq1qth5rX3zzTeL/dZbbxV7165d1bgcm0Nn1LoIrLmW5/a8884r9qWXXlrsJUuWVOO49tIfb7zxRjXu7bffLna+Fuh/nm9uJcz6G3k/xu/n/mnU+32facXbqLWncgwT+oP+zM80rdoWrHfDWlPZ10O898019ONXvvKVYv/whz+sxn3zm98sNmP9pZdeqsY9/fTTxd60aVN1jHHEa6ZVn7HVyntSaxi14pS+YQ3MiLp+3+WXX15stnCPqOsN5Wd+fmZ9mzz/rOH61FNPVcf4mTWQ8rrLOmMzqd1opo2IiIiIiIiISA/xpY2IiIiIiIiISA+ZE3lUlizx8wknnFBstteLiFi4cGGxcxtofmbLxNzSkpIopsWdccYZ1Timd+cUSaauMfWNkpuIiLvvvrvYDz74YHXstddeK3ZXm9OIOvVxXFLfWq3OeSynAFJ21kpRZSoc0/pzS0ymhefzYNox53+S0g0jRpc9zUab7FaLWq4BOWbpx9yy9JVXXil2VwppxMx+S6sl49DIMhim8v/0pz8tNuVQEd1ymVYctVocUraTJatca3NLTKbvD10Kl31FGQYlUBF1ejfvn1mWxNTdxx9/vNg5DZypuzmeec+kNPXUU0+txjHWW61N6dNxlQDkOfqI/HtmIu1s3ReZwk2/R9Tp4tk3TAOnHDEzyffFllSR83nJJZdUxyhVvOiii4qdZRKUJ7766qvFzpI07nfydzD+uKbmPRLjjxKcqcYOiZbMsBWLjFuuTy0ZA78/X1ddEvGIOha57+laUw7H0GO2RZ533qtuvvnmYjNGI+r9DdfD3/3ud9W4Rx99tNg5jjjvXJfz2kv/5/sDP4/jM+FM6ZIXR9TvByh7+slPflKNW716dbH5fiFfE60Y5t+mn/J3UOrK9wv5b3Ofw2fRfGwmmGkjIiIiIiIiItJDfGkjIiIiIiIiItJDfGkjIiIiIiIiItJDZq2mDbVpuQ4CtWm0c02DlStXFjvrxagPpJ43t//i36bGLGuTW7pW/hZqiVkvJ5/v9u3bq2PUJ7Mt9kxafPWZ6bT+7dL+tmqhUGuY55+actbDiKhbmLZq8Ewas92+c1T9NeMy12CgTjS3x6PWm3arpk2LSapjQ9+wDlhE3fqS9Rfy2s3vYBxlXS5rDz3yyCPVMdZN4fznv8XW77neAtvhcj0dSjxzXqiLj4i47LLLin3LLbdUx3gPZavtxx57rBr329/+ttg7duwodq4hRR/nGgysVcOaAWeeeWY1jnU08nXCegC8tlpx2Wcfd91bWvfFFl01x/Jn+ia3/D733HOLnVsQs803fZ9rD/V5zmeD/Pu6atrkWjKLFy8u9jXXXFMdW7NmTbG53rJuTUTEE088UWzGaR7He1yOsa49cF5TeSzHc9f1OnTft2oWcb/J9Yn1bSK69yJ5HeP1k31I33DO8152kvYsMyVf23w2Y923/Ly4devWYv/mN78p9tq1a6tx3I/kGkm8hmjnujW8FnL9lqHH3Efk5wWuV6zdFxFx4403Fpt1F1etWlWN476EsfjOO+9U41hLLMcUa97ynLKf+LfyvYH+5l5stu+tZtqIiIiIiIiIiPQQX9qIiIiIiIiIiPSQOWn53Wp71xrHVOqcVs10RB7bsmVLNY5pbExtyu3FTzzxxGLn1mxsL8704tySkSlQLZnIJKWeklYKMskpaEyTu/7664v9ve99rxrHdNPNmzdXxyjRGJokbTrMNE2/RZcf839nyirTudluOqKWvb388svVMaY4sh3xTOVRkwTTPK+44orq2FVXXVVspnxmutJN169fX4274447is2U44h6/aevswT27LPPLnZeT7nmU+IxFL8zVnIaPddApn1H1L578skni33XXXdV4yjfZRy1Wm3nlHOmllOSQzlORC2Boqwtor6ftiSy4yJpncv283k95X2Se5QLLrigGsdU7yx/455lz549xeY1ETG3v6sPtFp+c73iPjEiYtGiRcXO9zHOO9eohx9+uBr3xz/+sdi832VJ6Kj711ab21br60lpM9x6HskSB84z56QlI23FCtfnXA6C18uHH35YbEqNI9oy0qH5ajpwXtj2PiLi29/+drEZs/v27avG/f73vy82pYpcJyPabd15HmzjzjU6or7WKPHO3zk0n3aVu4io440y/YiIW2+9tdjcK+Y9C9faTZs2FfvXv/51NY4y/Ysvvrg69qMf/ajYvJ/SnxG137KM8bnnnis2pa75XUZrzzUKZtqIiIiIiIiIiPQQX9qIiIiIiIiIiPSQOZFH5dSv9957r9hMt8+dYphGlau1M82QVaCZkhRRp2aT3EWBkqictshuHcuWLSt2lke9+OKLU55TRD0HrVTvoaXCtaB/6Y8sk7jpppuK/Z3vfKfYuWo45zinFh84cKDYrXS0USvzT5KfZkqr0jpTRdklKKK+Ft58883qGOUVXBNySnLLP10dzUaV740L+fyZeso1LaLu/kNypXvGEbueUA4VEfH0008XO6eDMn2f10HuFsBjuSMOpVkbNmzoPN9xgv5irGTZxaWXXlrs3P2CqcH0z/PPP1+Na8VO1znl1GDeJy+88MJiZz9SapNlT13dWVpMkiSg1dmGMUtf5P0L0/D3799fHXvqqaeKTXlU3rONizxttuia9xwDlC5S8htR7225N1y3bl017rXXXis2ZTetbjM5xrh/5frKvXZEvR/Ox3i+Q9ijtjqXEso+OXcR3fKoPAdcxziPeU1bvnx5sfO6zjnn+p/lyqOuk5MG9428H0XUUhtKz3JXKEqiGB+tvWFelylV5FrMZ8eIei+V6epe1xo3LnHJ6zdLrrmOZZk111c+e+cSCvfdd1+x77zzzmJznY2o/ZZjrEtyn8sw8Pkkr+ss0cF1N+9RP6n02NVARERERERERKSH+NJGRERERERERKSH+NJGRERERERERKSHzFpNG+q0cj0aasR4rKUbzC3XqLnetWtXsbNmm+O6aqhE1Fri3E6ROja29XrmmWeqcaynk1ubdmmVx0WHOBu0Wi1Sy3j11VdX49iOmPVucs0U1rnIbRJ5nVFTmb+j1Y69S3s4ST48HJzbrPXlMdZjyG1UWQMl14YatT0xGbW+zXSOjQNZ+846CEuWLOkcyznONYUefPDBYv/hD38odm7rTd1urgPB+jRstZg132wHntcOxvezzz5b7FyLY5x8SB+wtgLrIETU85LXL/qL85Jb1HbNS55nar3ZKjUi4oYbbig264vlOnL821kT3vW3Wy2Yh95+mvB3Z/0/rxHq/0877bRqHNdJ1q2JqPcsrAP3SduQjjtd9VByvLGmVG4XzXV0x44dxc7xweuZ33/MMcdU4xYuXFhs1rWKqOOPe1vWFouo14e8JrDWxzitm6PQ2tdxn5LrErHeHv9dHpefTz6CbbwjIlavXl3s3AKaNVTYsjrfP7kOtGpPDZ18j+A+Mtfs432Msci6bxH1812rJlLrmmFsrly5stj53sd249mPHDuEWOyqy5OfETiXrfhgHcOHHnqoGvenP/2p2HwfkO9pXGtzHdWuuM/74fXr1xf7/vvvr46xpfioNQRngpk2IiIiIiIiIiI9xJc2IiIiIiIiIiI9ZE5afud0oJwK9hE5XXPv3r3FpvQlj2V7rlZ6PNPPc8oZU0op1Ymo0xMp18iSAMq0mGaXz2uSUhhbcM4p3TjllFOqcfQ1091yuiHnnymQGUpzsi+YQsd04fy5qz1mxDDSGWcKf3ueF6ZI0t9Zqsh2fllmmKWW0k2WUzCucio/07uZ1plTPh9++OFi79y5s9i5rTdT+3OKOFOXL7nkks5xXIdzaislVscff/yUvyNivGKR9yemCWcpG2Ux+fexjSjTifO1QP8wTTjHIqVYlENFRFx++eXF5j0yt+DkvTu3u2R6NM8xr+3j5Me5IssdeR9bsWJFsbOvOeeUQ0XU98x8vyNDn//W72vJ6rku5WubayLXpdxWmjFG8j6IbYyvvfba6hjlcdyj5n0oz2O2W8+OK10y/Yi6ZTPjLT9n8DPliZRZ5O/I80853eLFi4u9ffv2ahz3wHk/NEmlF7K0hvOepTXcU27btq3Yue0245v7ivy3FixYUGzuYSJqqSLvY9xXRdQyybx/YiwOTcLfWk+5L8mSKrbs3rhx45R2RD2X9FuWtlLufeutt1bHuA7w2slyU5YI2LRpU3WM8TyXfjLTRkRERERERESkh/jSRkRERERERESkh8yJPCrD1C+m8+UUMabrMkUpH2tJj5hixfTiXJn/zDPPLHbu1sGxlEQ9//zz1Tim2uXUx67q9eOY3jZTcrobU+OY0j1qGmFOZWXaY0sexcr8WSbC6ydL8viZ0oPp+HCI/u66nnMsMv5YrZ2V/SO65XARw+5wMRu0Op0wNZvXb0Qdc5RA3XvvvdU4+oNrXCu2c4o4U4uZgpzXZK4J+Vrisa7OBOMG44Ndaei3iPr3tmS+jDFKiCPq1F3OZU4rv+iii4p98803V8e4jvK6yJJGdiXK58v1l+fe6qYxqeS5oySNvs4yKvqa98iI+p42SdKKw9H1+7Mchdc60/cj6jWL13nukscUfnbAybHIf5clGVwvKE/M8qhWF5OhrKOfhLzWUkJBSW6Wy3AN5bqY72lcC/P+kvsgxnPubMPOifk+3tV1aCjx3OqmR0lxPsbnBs5Ffobg91OeSJlcRC2Byt2G2WGI621eH/hMO85dL6dLS7LENSnvWfh8TQlojkX6lFLU3OXr+9//frHZoTifF58l83748ccfL3Zea4/U/dRMGxERERERERGRHuJLGxERERERERGRHuJLGxERERERERGRHnJEatpQ39XSfbF+RdbbzkQvxu+gdjgi4qtf/Wqxs16YWjrqhXO9jZZeWNo6YPo6az9Z04Y+zLVQNm/eXOzsG9ZLoFaVmuWIWoOcz4O+b9Viou/zddtqqzoEWjpq6sXPPvvsYudYpNY3z21X3A9ZAzwdeL1lPT1rYOSaW9TJr1u3rtijtlzPdTSo82YNm4g6FunPHCv0aa5pwutiXGsxtM6bvzfXc2MthKynZ4ytWbOm2FnPzbpR/I6FCxdW41jHIbcm5jW0d+/eYlPTH3FoDBNqx2m3roVJovW7GVf0YZ67ffv2FTvf03g/GrXN7CTQtfbkWGQL9RyLrGNDn+QaJdzHtGowsDYYa9hEdNeCyHVTGItceyeJfG2zngjXxQz9kWvf0Ifc1+aaM1wb8z2Y/4737lx3hb7Pe2B+56jPIOMU662aNqyPl1tJ8zP3/Hkc5/PSSy8t9pIlS6px/Nu7d++uju3cubPY3FfxOSaivjZyTZtJeX7MaxDnhHVrIuo5op1jkZ/ppyuuuKIa961vfavY8+fPr46x9tSLL75Y7A0bNlTj+JyZn+2OVFyZaSMiIiIiIiIi0kN8aSMiIiIiIiIi0kOOiDyKtFJyW+njM5FEMVVq9erV1bhbb7212DnNla29mQqXUylNLz7UZ5RNUJ4REXHaaacVm2mKWQrBtLM9e/YUO7dYY+p3/g7+bbZ3z9dBV+vMiFrexe/PcgBeBzntcYjQ5/R3Tt1dvnx5sc8999xiZxkP07vz/I0qjzL+Dm1nn1vIEsqemFI/aqvlnKpMmU2WAzAWeY3kNFd+Z04l5zrQkt/0mXyN8tpmDGzatKkaxzT97FPKNxiLXPMi6jWKa16+93HNe/fdd6tjXH/ZMj6vh62U/VFb1E5qPNOHOT7OP//8YjOO8prJPQvTviO65VGTTpc8KsuNuDfM8cG44j0uy4EZY1yzKSGOqFtJ046oY47p/JQtRtTt3/PaPin+z2sJ732cu4haKsxYzDJS+pAx9sILL1TjKJHJazJl+zlOCe+LueQAz5F7gSGun1mSTbIsjfdJ3uOyBJESYEqi8trLZ7+8/3j11VeLzbjM47j2znSfOy50lUTJcnuWImnJNzlf+ZmT8UEJ8Q033FCNo9w071/5LMn9F30bMTM54mxjpo2IiIiIiIiISA/xpY2IiIiIiIiISA/xpY2IiIiIiIiISA854jVtWno9atVGbQGax/Ez9Yo/+MEPqnFLly4tdm41Rp0raw1kTfC46g0/KaO24cs64BUrVhSbet6sDacOm5pgtjKNqDWjuZ4KNa2sY8O2fhG1pjzrXVkvgtrLrEcdelvNVt0iavdzGz223GOdkxxH1OG36kbNNrNRQ6tPZL07YzHXEeJ1z7azWTfOOeL357beV155ZbEZ5xF1C1Pq+nMNHsZR1hI/8cQTxaZ+fRz99BFcK6iFX7t2bTWOdTVyTZuuuls5xviZ10JeN1l3LNckY90FrsWttrn5PDiW+vAhrpuj0lUXLN8/u2ra5HoYrP+Ua7J01awa5ziabThH+drmdZ/3LfQJY4yt7SPqueYamFsJc13Oew7WYGCdndYeKbeo7arJOIRrgb8h7yFYV4OxEhHxyCOPFJt7Pq6LEfXekN9Bv0TUe6I1a9ZUx/h8wnPK8Uzf598y6nPRuPq0q9ZURP1stmvXruoYa9Jw/5GfV/idjJXcGpzPIbn2JddY3o+zr+jHViwOjVFr2vC5L8PrIF/brFnE2pkXXHBBNY5rbT4P7jf5/M84j+iHn8y0ERERERERERHpIb60ERERERERERHpIUdcHtViNlL4mNJ9yy23FPuyyy6rxjElmemlERHbtm0rNlNg+5Aa1QeYAsg2lxF120q2fI6IWLlyZbEpk8jtfXfs2FFsSgVySiHTkXOb4auuuqrYlERlWQdTl3MKMn8b/1ZOneR5ZYnKuNLV1juiTjFlCveiRYuqceedd96U/yZLX5henNMWGXNd6dzToSWJGke65ieinvOcFsxUUcbOG2+8UY2jP9h29vbbb6/G3XbbbcXOMkOuyYy3nA7L1P577rmnOsaUZMbpuKZ9R9T+ouSL7bQj6vRrpnpH1HObJXCEayzTwLNEjddC/luMHcpGsrSJPsltT3ke9OMkSY/zesqY4Hp6zjnnVOMojeOc55jduXNnsXN6N+d5kiVpLVqSDF7PLekffdxqVUzf533Q1q1bi51lbpRTbt++vdhvv/12Na7VonbIMUby76TfshSJscM1Oa+tnGdKZ/Lf4uccp9wr0k/0bUS7bXuXD4fiW/6OvD/nPNFvEd3SwrwP4ndwPcxSf+5N6O+IOoYZ6/l86eOWH4fiu4+Y6XraRX7+ojzq2muvLXaWNJK8nm7ZsqXYfB+QfdgHzLQREREREREREekhvrQREREREREREekhvZJHzYRcmZ8dF6677rpi5/RGVht/9NFHq2NMxWdK29DS1mYK0+TZBSqi7nBx0UUXVceWLVtW7LPOOqvYeV75nfRF7uRAuQY7REXU3WyYJpevF6ZE5lRZppYzVTafxxC6oLQ6ROU5o1SMPmBHtog6nZ9z+dxzz1XjKJfKKeKcTyVRh9JKH+7q8hVRS1+uueaaYi9ZsqQaR6nT1VdfXexLLrmkGsdOVVkOwHOk3CrL5O68885iP/DAA9UxxuYQZapMDc6yMf7eLB/kvGepE+G/4/e3/k2Wo7IjB32c05W5Po6aIj70eyvnK0to6QPeq7I8iusu5RnsdhFRy4uzPKpLWjj0+Z8peV5a93fGaasTKj8zBnJssxNRPsb1kNdC6/7ZkkcN2f/5t3Ee8vrUdZ/JPqQ/aOf7LGM9r6f8W7RzzI4qj2r5cAj+zbHHPXmWnnFvQUlU3stSRsXvyLHCTlWZLh9naTCPDXEPMxM4D/na7upcmvcs7FLLjlH5Pst9SX4Guf/++4vN66CPcWOmjYiIiIiIiIhID/GljYiIiIiIiIhID/GljYiIiIiIiIhID+ltTZusJetqQUxNf0TdVpp6xdwS7i9/+UuxN27cWB1jO7BxrVEyl7TqaFBHeMYZZ1TH2E6WtTKy9pDaxq9//evFvvDCC6tx9G9u+c123Tzf3K6PrfyefPLJ6hhbbnJcrjlBTfNQtKqMt9wmkfNOH7Pde0R9bbA2UY63119/vdhZuz9qTZvZqFvTR/3qdMg1mTZv3lzsFStWVMfOPPPMYrPWFGsURdT1pej3HLOc/zyPrLPAOg2/+tWvqnF33XVXsXNLxnH3zeFotcWkFjuvL4wxauhzXSHq6RlT2Y+8hriGRtTXwoknnjjl342ofZfrf/F8+VtarXLHkVYdk7yeslYQ69iw7ltEdz2oTZs2VeNadeBmo0bYJMM5a93rWzWMumqNsWZRRL1u5vviqPsMfXwonLfWnh0AAAkPSURBVJO8v++qYZl9yPlnnZT8PML6G3lN7nqmyefUtXbncxzlv48z+Zrn2sb9eUS9VyGs5RdR+5UxxjpREfU9LdccYttwPl9MZy87RH9NRau+VL5n0jesm5drvbGmDcflPSTr2Nxxxx3VsfXr1xe7j22+iZk2IiIiIiIiIiI9xJc2IiIiIiIiIiI9pLfyqBZML2aaf0SdUsxUb6blR9TSgdx6NrdNlBqm+WW5ESVFixcvro7RV0wpzWmP/E6mD8+fP7/znHIa+CuvvFJstkTdu3dvNY7HXnjhherYnj17is1rKbduHGLKeVe7vYg69bQlc+P8tdrtMY0xz21LQvFJGYqvPiKn4z722GPFXrhwYXXsu9/9brEXLFhQbPozok7Zz+0ySUsyyRbEv/jFL4p9zz33VOPYwnNovjkcrd/L9SVLkbpaBrfkay0/jtomntKpfA9oSeXmMp77TGte2eabcZpT/DnPbEG7ffv2aty7775b7CO5nk4CrRbL/Ex/M2U/opbMMKayvJjXAiV0EfX9lNdCju3ZkA0Pmbz35HpKf+Z55frKe2RLUnryySd3HuM9s3VdTefY0OGcHThwoHMc18NcsoH3U8ZplipSfrVt27bqGOWolFXl50hLbRxKS27K53zG1apVq6pxXBv3799f7NwG/s477yz2vffeWx2j5K3vcWSmjYiIiIiIiIhID/GljYiIiIiIiIhID/GljYiIiIiIiIhID+ltTZusxaWmlFrQpUuXVuNYN4WtmXfv3l2Now4xt8pVe9iGmj/WoYiIeOaZZ4qd6whRR0j9dm6FyO9stV9r6SGp+aamNZ8v66nkmiCsBzDJtQCyf6j9ZVzmFnucd7b1zrWDqAPO/p6NeZ8Uf+UYoM47tzhkzacf//jHxb744ourcdThsxZHnlNqgu+///7q2C9/+ctis56RtcOmZtT6NhG1z7le5Zo2jNPWuslYz3U0qCun3jz/Lfo1ny/v65NUb4PzmtsCsw7fokWLin3cccdV43hP434m32dbNW0mZS2cK1rXLH3M+MhxxLoay5YtK/Z5551XjeNelrVWIup6cYxta9p8Mrragbfihn5nvaKIumZRrhfHa4R7z7xmdtXZOdx5DY1WfbRc6421K7ke5vpr3HuydXfe87I+Sq6LyX0vnzmtJ3Z4WrUzef/j82Ku/UVYo/SBBx6ojrGGIuuARRy6D+ozZtqIiIiIiIiIiPQQX9qIiIiIiIiIiPSQXsmjmCrF1MGIOrWQKfzXXXddNY4pppTC5PQ5HlMONXNyWhnTPHMaIdvmrVu3bqTv5zWRUxb5uZUGzHPM52vK4v/TSj3NaZ6UwlBmk9uud10LlEpF1CmlOQ18kltaflLowyxdu/vuu4v90EMPFXv+/PnVuOXLlxebaf1ZSrh169Ziv/TSS9UxXhfjlIY6DnSl8+d55jimIec1len9+R7Meyivp5Y8p9X2dJLW3tY6xvWVdk7lZ/r+pk2bik2pd0S9nua9jevp9Mj7ilY6P6VJjKO8pl511VXFXrNmTbHZ4juijo9du3ZVx+hjxmW+f8r0GDU+up4ZKCGOqKXkeS2kD1miIT+r8G+5P/qYlsyXMnvOX55byqN27txZbPotf1/+Du6F6OPWPXhSyesp9xiU4kdEnHvuucW+8sori005cUQ95xs3biz22rVrq3G8n47zPtRMGxERERERERGRHuJLGxERERERERGRHnJU5VE5VerYY48tNlPxIyIuv/zyYt92223FXrlyZTWOaamUZOS/1UqPmuSUw7lkJnPZSoGUuYNznaUw+/fvLzalEFlqwVReppfmFN9Wlwbjb3bI88g5Z2p27qS3Y8eOKb8vr6f6ae6YjbmlvynBybFN6SM71ETU91NK3rL0juNyt76ujnxDoxVvTMmPiNiyZUux2aUkdwKibyiVynPM9XXIc3w04D0uy6OY6s99aO4axI5C7CyV117eZ3PHxWeffbbYr776arFb14Jr9Mxpyccpz8hdaSi5oRwqor6HMrazjIp7JyU3U9OaB8ZAlrVxrumDHNsk+6AlS5aa/IzATornnHNOdYzS0csuu6zY7B4dUccc74t5PR1KrJhpIyIiIiIiIiLSQ3xpIyIiIiIiIiLSQ3xpIyIiIiIiIiLSQ3pV04Y6wqxbW7RoUbEXLFhQ7JNOOqkax1Zt1MtlLWNXvY0IdYkiXa2EM6O2U5+k9r6TgD7sN6PWYMhtpVnjJrfy5v2a90/Wqcmf87GudWDo1xPX0Fyzgpp81ghrtQZv1WmQ2aNVm6gVYzyW69HwGGMsxwqvi+eff7461lVfyhbvc0PrOmALaNakyuNYbyOiroHE78hrMo+1/CuHp1UzkfGb6y7K7JDnn8/e+X0Aa7pxjct1+LZv317srVu3FvvAgQPVOK6v4xw3ZtqIiIiIiIiIiPQQX9qIiIiIiIiIiPSQoyqPaqVK5XZdTBlki9Hc4pAtxdgWkWlTEXV6aU5LFZl0mKqY2/SRmaSLi8iRpUsmkaXBTAtvtXUfNZ6N+0PJMlL6oJXCrdTl6NOSDfMz/Uh5S0Sdtr9+/fopvzt/R45TJXFHl641NMtq+HyS91FcX1vXlXEvQyHf+95///1ib968uTpG6ehDDz3U+Z379u0rNtfWoUrczLQREREREREREekhvrQREREREREREekhvrQREREREREREekh86ajkZw3b96cCiqp+cx6+mOOOabYbAd+7LHHVuOoIaWd22y2anH0kYMHD847/KjDM9c+lCYbDx48uHo2vkg/Hj2MxUFgLA4AY3EQGIsDwFgcBMbiADAWB8GUsWimjYiIiIiIiIhID/GljYiIiIiIiIhID5luy+8DEbF7Lk4k4tB2YOSDDz6Y0p4Qzp7F75pTH0oT/Tj+6MNhoB/HH304DPTj+KMPh4F+HH/04TCY0o/TqmkjIiIiIiIiIiJHBuVRIiIiIiIiIiI9xJc2IiIiIiIiIiI9xJc2IiIiIiIiIiI9xJc2IiIiIiIiIiI9xJc2IiIiIiIiIiI9xJc2IiIiIiIiIiI9xJc2IiIiIiIiIiI9xJc2IiIiIiIiIiI9xJc2IiIiIiIiIiI95P8A4Sk/UtxFrqcAAAAASUVORK5CYII=\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: 13.7 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,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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0244, PSNR (dB) = 16.12\n",
      "MSE = 0.0586, PSNR (dB) = 12.32\n",
      "MSE = 0.0114, PSNR (dB) = 19.42\n",
      "MSE = 0.0214, PSNR (dB) = 16.71\n",
      "MSE = 0.0331, PSNR (dB) = 14.80\n",
      "MSE = 0.0110, PSNR (dB) = 19.60\n",
      "MSE = 0.0476, PSNR (dB) = 13.23\n",
      "MSE = 0.0471, PSNR (dB) = 13.27\n",
      "MSE = 0.0768, PSNR (dB) = 11.15\n",
      "MSE = 0.0338, PSNR (dB) = 14.71\n",
      "Wall time: 3 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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5447\n",
      "0.3796\n",
      "0.5834\n",
      "0.7767\n",
      "0.4149\n",
      "0.6362\n",
      "0.2455\n",
      "0.3286\n",
      "0.2087\n",
      "0.6214\n",
      "Wall time: 188 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": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8544\n",
      "0.6382\n",
      "0.9411\n",
      "0.9042\n",
      "0.6124\n",
      "0.9616\n",
      "0.6467\n",
      "0.5846\n",
      "0.5175\n",
      "0.7892\n",
      "Wall time: 158 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": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Simple Autoencoders in Keras - 784-8-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
}