Simple Autoencoders in Keras - 784-1-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 = 1"
   ]
  },
  {
   "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": [
      "1595258315\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 22:28 - loss: 0.6930WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (6.180007). Check your callbacks.\n",
      "  3/235 [..............................] - ETA: 15:00 - loss: 0.6928WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.766011). Check your callbacks.\n",
      "WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.384506). Check your callbacks.\n",
      "235/235 [==============================] - 46s 198ms/step - loss: 0.4261 - val_loss: 0.2889\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 23s 97ms/step - loss: 0.2804 - val_loss: 0.2737\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 26s 112ms/step - loss: 0.2712 - val_loss: 0.2680\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 26s 111ms/step - loss: 0.2644 - val_loss: 0.2631\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 25s 109ms/step - loss: 0.2626 - val_loss: 0.2616\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 19s 82ms/step - loss: 0.2449 - val_loss: 0.2441\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 18s 77ms/step - loss: 0.2445 - val_loss: 0.2437\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 19s 81ms/step - loss: 0.2441 - val_loss: 0.2434\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 18s 76ms/step - loss: 0.2438 - val_loss: 0.2431\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 18s 76ms/step - loss: 0.2436 - val_loss: 0.2428\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 17s 71ms/step - loss: 0.2434 - val_loss: 0.2427\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 16s 69ms/step - loss: 0.2432 - val_loss: 0.2425\n",
      "1595258315\n",
      "Wall time: 19min 23s\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": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAADiElEQVR4nO3dQWrDMBBAUankCl33/sfKPndQD9DQEnDjb/He2gjBMJuPjedaawAAAADQ8nH2BQAAAAD4SbQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACLq98vCc0//BT7LWmkecY4aneqy1Po84yBzPYxe3YBc3YBe3YBc3YBe3YBc3YBe38HQXvWkD73M/+wLAGMMuQoVdhAa7CA1Pd1G0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIEi0AQAAAAgSbQAAAACCRBsAAACAINEGAAAAIOj24vOPMcb9Py7Cr74OPMsMz2OO12eGezDH6zPDPZjj9ZnhHszx+sxwD0/nONda774IAAAAAH/weRQAAABAkGgDAAAAECTaAAAAAASJNgAAAABBog0AAABAkGgDAAAAECTaAAAAAASJNgAAAABBog0AAABA0DeT5kt68dYzGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2daZfd1NFGi8xzgBiMzWACJARI/v+/CJkWWYt58MDkhCSEzPH74V0o+zz0Kasvdrd0e+9P51pqXV2Vqs6RXE/VA3fu3CkREREREREREdkWXznvExARERERERERkS/iSxsRERERERERkQ3iSxsRERERERERkQ3iSxsRERERERERkQ3iSxsRERERERERkQ3iSxsRERERERERkQ3ytdPs/MADD9gf/Jy4c+fOA/fiONrwXPn4zp07j9yLA2nH80NfPAr0xSNAXzwK9MUjQF88CvTFI0BfPApO9EUzbUTOjnfO+wREpKr0RZGtoC+KbAN9UWQbnOiLp8q0ERER6XjggfE/ee7c8T9rREREREQOxUwbEREREREREZEN4ksbEREREREREZEN4ksbEREREREREZENcmFr2mTdhTVYm+HLwWvO8Ve+Mr475Oduv86GtNV///vfE8f5Oe3Lz9r+bFjrl9rjyzHzq7z+3Nbtdy98cbZfbtP2JzOLr1WHxdSuNlEXN9faWOQiMouVp4mpM1zD3B8623Qx85D1TBdPteHpWWufQ581yNp5sfNTuTccEj/vxhbsZKaNiIiIiIiIiMgG8aWNiIiIiIiIiMgG2b08qkuBOiQtLunSFrt0/otKl2L41a9+dRl//etfP3Gcn7/xjW9M9+Px0jb/+te/lvE///nPZfyPf/xj2I/b/v3vfw/b/vOf/yxjU1RPz9oU4tnf5GdTiE/HTNpUNfrO1772v2mA/lY198Xcj8fvfJH+R9/LbZ0vXuS420lJaUeOq+bxNu3I+yKhTWbxNT9zv6q5HfXZ03FoPJ2hrObLc4jUgv6WvtfJUWcxkP5VNfqsktO7M5sz186fGXe7eMprTrtlzKQNc1503fP/dHLgzsdor+6ZpLMjbUDbpR27Z41OUiwj90KqeOgz/xYkbmbaiIiIiIiIiIhsEF/aiIiIiIiIiIhsEF/aiIiIiIiIiIhskN3UtJlphFNrOtPrf+tb35rul1q0mb70NBrFY9bud3rtrj7Gd7/73WX8gx/8YBk/9NBDw34PP/zwiWP+fdVow6yr8Oc//3kZ/+EPfzhxXFV1+/btZfyXv/xl2Pa3v/1tGdP2p2kbfozMtKFr9cLf/OY3h/1SP0xmeu70RX5OXb+++MVrTl/64Q9/uIwffPDBYb8f/ehHy5h+mr7I7zpLX0xt8jH44tpaYFVV3/72t5fx9773vWWcdmQcfeSRR5YxbV81xuz0sZkdabf8vNaO6bMXVde/tiZcN89mXQ3C69rVyqA9OttcNDsdUrOvi73f//73lzH9t+qLdiX0o7/+9a/L+LPPPhv247a///3vwzb6H+1/7O2IZ/VOqka78Zkh57vZ+jXjKY+RfkTb/OlPf1rGjLO5Le3LOnDdGvUYbNjVrcl5kb4z87eqcV7kWifnTx4v7cg57o9//OOJ4/z86aefDtsOedY4NroaNIf4bD7zd/Mibcr166H1+u7nOtRMGxERERERERGRDeJLGxERERERERGRDbIpedS9SA3+zne+s4yZ4pYSnExFJUwjZRob0xRzW6aeztKo9poGPpPBrL3+VVWXL19exk888cQyfvrpp4f9rl27duLfMCW1arwP8vozff/69evL+O233x724+f33ntv2DZL88+UuWOT3FT17fFmLRTzXqB0g2mpTEmt+mLqMWH67yyduGpMG16bBn4MvjiToFWN1zXj36OPPrqMn3zyyWWcvvjUU0+d+DcZP3kfMNW3quqTTz5Zxjdu3FjGnS/SZ6uqPv7442XMuMv7o6r3vy37ZieBYppvxsBLly4t48cff3wZpx2fffbZZcz4yr+vGn027ch4+M477yzj119/fdjvzTffXMZpx48++mgZM6amHbtWxXunS/OftXCvGmU29L+8JzgH5zqKcxdjJmNr1Wj7nO9msqpjTOXv1qE539FPO6niY489toyvXr164r9XjXNmzlWU0HzwwQfL+NatW8N+H3744TJOuQb9j3Nml/a/5RhKOtk24yvjXdVoN65T0jacFxlrr1y5Mj1eShDXrlG5LqU9q8a5lT7c+exebFi1Ph7mGpI+R9txrVNV9cwzzyzjH//4x8s47ciYmteWaxPOi2+88caw31tvvbWMb968OWzj3Nqtb45tXuxKnczWRCl7mvks5W65X8K5kM8W9K+q0TZrpYpduYZDMNNGRERERERERGSD+NJGRERERERERGSDnKs8aq0EqmqeHpXVwJkSxVS4TIvjfplmxjTS999/fxl3KftZ8X0mL+m63mwpbXFtpXamDVaN6WlM/a0a0w+ff/75ZfzTn/502I+2oqwjZQNMM0tJzKwSfN4vmWpHeF+slbjttfNC+mKXljrrSpQpqrwXKIfrfDGvF1OImRpMmU3ul77IeyPva7JVmVsnT6MtTuOLTAum/6Uv8u/YGSNTzjtf5HlQypH3C30xf/NMkpEp53vqHjXzsUzZ53XPNH3G1J/97GfL+KWXXhr2+8lPfrKMKaNKaQ3PqZO58TxS/tHJc3ifrJUq7jGmrpVAVc0ljjk30VZdyj+7g+W1og0ZT7nOyf1SOsU0cP7O9MWtxtO70XWl6eY7+illh2kfymmee+65Zcw5smpcq6Qkg/Md5TMpVexi6ozO375sav/9ZCaJyu5djE9rZcOMn/mZ405W00n4KavJc6INc96ddXLrpIpbtmHVfF7s4mHKfGk7SoNffPHFYT/OmZxL83i8h3JepB3pw2lH3gudHdd27tvLvHhIx72q8ZpzTZQd2uizlC1yXDWuQ/O60obdMz+lqPlbZs/8XTmNQyRuZtqIiIiIiIiIiGwQX9qIiIiIiIiIiGwQX9qIiIiIiIiIiGyQM6lpM9O0ra1bUzVvodhp/Kk15b9XjXrD1LexLSnbl2atAf6W1HPP9Pr5XTzGlrSmXR2Ntbp7auurRr0n6yrkfjwm9aPZspLb8vrP2iKnzpT3T7Z34/exDVxqWrv22FvVmVbNW3dXjTZOTTi1/LQ3taVVo3afNYyoMa7qtabUkHa+yG2dL+5R9722vlT6ImsipEabtWqow0//4PFZ2yLbHfJzXv9ZPY+shXL58uVl3Pkiaxat1RVXnb8vdjGV7YOzVgZrPmXdC9bEoD6f/1412pj+nC1FWbsr7cjzpe3ynNgWM2P2zI5dTN0LvBczns7q1lTNW0Wnz1KjT1tnHSr+XRdPX3vttWX86quvDvuxPW3XsrTT59OG5+17d+OQmMp5q2r0MdoqfZHzH2tvpL353WkDfjfnwoyH9OGM2fQ5xoGuXfSW4W9nPM1ab5zjsgYNnxPoVy+88MKwH23INWTWCGMcyOs4qzeUtmYtHMbWqrFt+8ye+d3pp+ftm936ppsXZzUTq0af+/nPf76Ms9bbtWvXTjwev7dqvEYZv2lz3k9Zw4j+161vaOO9zouz+lJd2/Zc08/mwvRZ1mek3WnbqrEWTvoA69Ny7sv4z/sg4yQ/d8/8tvwWERERERERETlCfGkjIiIiIiIiIrJBzlweNUsJrxrTo5gaVTVvocjUqKpRkkEJDlPvq8aUya5VJVPaMtX71q1byzjTvmbH2zKdtGCW+t1JZ1L+wM+8/tkGna1IKUtiClvVmEaYKYtMY+Z35e/itkxt5TF4/K3JLk7DzI7pi0wLXCvXSAkifZM+m5JGHj/vBdJJ5dimL1PEySEt9s6Drk3iWhvyumabRN7rtHWmVfO6MqWXEtKq9b6YcZ0w/ne+yN+59XThQ9pdnkZyylRhyhMz1Zjz2O3bt5cxJUpVowQu5zQek9vyvqOt8r7jPcljbN2OM2a+mD7QyTVm6fW5tmHqN2NttpSmj2U85TVnPM10/a616dqWpVu26aHyb15bzoNV41xIaX7akWtRHj9txXVpt8bg/ZRrLvpi/pZuzboH1sbTjIWdpJuyNkqgunbsnDNzXqQNO3+gLXjsPN9cizGudOueLZPnPZOSdnbM57tZ6+f0WdqE82JKm2jjPF+eI6UvXZzP9c1M4rjlGNqxdo3atfLmWofP9RlPKYPis0XaupOb0k/5vWwFXlV18+bNZZzxk7a6n88Z+/RyEREREREREZEjx5c2IiIiIiIiIiIbxJc2IiIiIiIiIiIb5L6IWlOHN9P1d7UaUn9LvRt1qKlbo86M7fDynGbt9qq+WNdhBnVxnQaSWuXU0m2pHkpX02Zmw9TudzViqMflMWinqrGODesGUU9YNdZpSK0v7xHWQ8r7qrvn+Hmt/ntL9qzq7bi2bkrWLZr5YrYs5b3BmidZm4h1NBL6Vaf1pSY1W/HNfLFr+X3edFrmmd1SL0xfTJ186q0/J2ucsNbFWl/MY/MeoV44NeqdL87qqXRs3Rdnc2FXm6irE0Y/zVah1OuzXglrFuXfZf0h2o7f27V6ztbXMx/eWhvaQ1hbEy7nTMZTavKz3gZ9gnExdfeMf3kdua2rKcR1StfalPE0v2tL8TRZu77pWtRmTOValOOMc/QxxtuuFWzGBJ5Ht76kDbJmztp5cau+uDaedmvUrKNBuzHGZRyjz3Ftk/NnV+OE9w/vs67GZtqG+65t670Fe6591lhbsy/nRdqR1zn9g/UzOUdyXDX6Std6fG19xmRmn7364iyeph/xvs8130MPPbSMufbgv+d3cb2a158xM593ujmT8P5ZG087fz4EM21ERERERERERDaIL21ERERERERERDbImff8Y7pXpg0x9SjTwphGxRQrpuVXjamK3I8pbFVjilWm3TF1iueYsinKOvI8mL7cpRBvNd2to2v5zVTErq0d7Zup/Ew9Zfo+0/qr+hQ0flfX8pvn0UmEunaKW7Zhd27cdhpfpP07X+RnpiinjIqprF0bPZ5jyhHp9ym34r78Xcfmi12r6EwDn/liXle2VmcqcfoibZMpyEwZ5jhTZXn9cxs/d7KOPdqwqpdk0FYpWeI2/vaUnNJ277777jJOO/JeyDRkxnqeY5dqnCnEs21bT+dfA88zfaCLOzM5QPri9evXlzElphnvmEqefsTzOHRtw307CdRe7JasbWme15brSNqUc1N+5ji/i/E8/Z5yEJ5v3jO0Y0p3+N2z9epJ57V31kriacNONsw23+mLXPekhId24/fmepj+l3Gdn2nDXM8dgw1znqG/pfRsNi+mHWk7Sr5Tcsp7IefFjAOfk3akfbq1MuNrJ3fccnydxdDuOSNZ+8yfNv2cXPPS//KZn9/F+NfNi+mLtHcnt/qyvmimjYiIiIiIiIjIBvGljYiIiIiIiIjIBrkv8qhOdsB0qK4DTAdTj5jeVjXvYpIVv0mmt/G8mB7FrkZVY4pkpkV2qd9bZW26HdMUszMCJRmdPIrXONNxmUZIW2cqHVPaUqbFdFPat+vK0923XVrinlL5196LnS/yGtI+2RWKvkjZXHbd4Dl1khn6YqZE8nOmT/K3bN0+p6WTRzE9tJNHzeSgVeM1p607yVymnvJzJ6vpug7NUpCTLfvi2jTZTh61tgtY+gDnLsoiMvbSVhlT+Zn26Loo5LzIz53MZo9zZieB4vXK30Y/6NY2nBeZ6p0dNGfHzu/m9edapmq8X1Jys5cU/Y7OF2nHTuaWfjrrAJRSC86T9J30ZX7u0vk7f6NdKXXNfbvYvlUbdzac2bNqnGeysxRjHP8u/YOy/Zk8I4/XzYv8LSmroS/mvcTv7uRRW7VhVW/H7rzpA+k7XAt18ij6Iq9zdszj8btnjdnauGqcd7t4S/vvxY5rn/nz96xdm3fPGbyutEXaiWTsnq1n8n5hDE3Za/esei8x00ZEREREREREZIP40kZEREREREREZIP40kZEREREREREZIOcScvvWVvlQ3XsXTstbnv00UeX8ZUrV4b9+DnrJcz0/6mlo74tW4PtRZM/o9Mo8nplTRtq7bPF4awNX/fd/K5se0n94uXLl4dttG/WgSBZj4F0ba9n57t1ZnV6ulbo3THob1n7gHpV2u6xxx4b9nv88cen50HtL/WlWe+B+6WWeE/2+Zy1tVk6/2BryvRFarS7OmOEsTu/i/r8zhdZWyfrbZCutlHH1mzd2XF2rqm3Zj22rk4Yr2dXi4O2evjhh4f9ePz0U9qVtSCyVgb9L2tszOrp7HW+nNmwW9t0NRxmNaSqRhtyTsu1De2U38W4SbuxJXzVGGv3UlfhNKytaZO/nT7W1UXgMbJGCY/JONz5W8ZbfvcHH3ywjLMNLdesWZ9hVlNqL7641ob5exi7shYK62DM6ulVjbGL69CsHffII48s46w9xe9mXMzvop9mrGXtsq4+yNZ8drYOvds2wuuX/pHPJZ+T633eJ7R9Ho+2y/UN11ldLSv6Zta0WTsvbs2OM2Y27FqY59qTsZbXJJ8zOE/STt0zf9ay4vWnT3XP/F2tt/uJmTYiIiIiIiIiIhvElzYiIiIiIiIiIhvkTORRZG26eJcG1rUIZrr3Cy+8cOK4qurSpUvLuGvlzXSolGQw9XsvKaUda68/U3M7SUam8jONlKllmTLHtDimQGZLY343JTZVYyocUyWzFS7PI9Odma7Hc+paFe8lfbFqvS92KZq0XaYcMlX4xRdfPHFcNaY0Ztro7du3lzHTTTOdn3bdkw1mHOKLKQOkL+Y2pgLz2nWtCulHaWv64tWrV4dt9E1KfVL+QbI9KuM8/S/j/5Z9cW3Kev52xtHOjl0aMvejv2VbTG5LuQbvJ6YrpwSK55Hp6LMU/k6Wtxe633OIHCDvA17/X/ziF8s41zZcAzF+VlXduHFjGTOG5tqGEo2t+dH9YGaD/O1r5VFd2j/9+YknnljGTz/99LAf5f255qDU6datW8s417L83En492jjtWuWvP6cuzL+0abdXEgbUnqcMZM2zLUyj8+5MOV0tHWuX3mMPdrwJGb3ZfoAnwdS5kYbd2sp/h3Xqykn72RuPA8+L+Z30a6dZPIY5sXZvXiaZ/7ZeqZ75n/++eeXcT5n8Jk/23XzM2Nm114848NZ+Z+ZNiIiIiIiIiIiG8SXNiIiIiIiIiIiG2RT8iimQ2XVfqZZ8++Ybl81ppsyPerZZ58d9mOqYqaNMlVx1r2mqk+f3DtdqhdTDzPlk2mjmbLItDbaN23N/WjfTIujVCCr9vO8eB6Z0kx5VNdJgH+315TFjkN9kaQvPvnkk8uYvvjMM88M+9FWmTZK32TqaXbJOMY04c/pfg9Tc1OOSJ9I28w6naRtuR/9I32ANszUYm6jz+Z9xXsuz5e/k/Gnkypujc6O/B2Zsk9ZWkpmSF5PMksDz/jNVOOu4xj9Ms+X9kkZHc+f99aW7XYIa1O9q+bxNKVw165dW8YvvfTSMk5ZDX0n4yTjK1O/L1I8/TJ0kiL6H/dL/6Bdn3vuuWWcdqQvZtegWbehlM90nVbpc53/7cX+s/PM39b97lnXqVw30jZPPfXUMk6ZPqU0GZ/pczz3TlKacN7Yawxde3+lDTi35LaZhD+vEedCyk8ph6oa7ZhzML9r7XNCrlv4d11XzS3Lv9dwmnlxdt/ncxpLYVAqnM/8XB+//fbbwzZ+N5/zU27a+eJZYaaNiIiIiIiIiMgG8aWNiIiIiIiIiMgG8aWNiIiIiIiIiMgGOfOaNqTTBLPWSNWo4eW21N1Tt8ZxavepIcyaNtSxsc1wntMx17RJqKWkRjtrT1BvmPUMaO9Zi7uq0Tb8rtyva5PIbTzHPF/qy1MrSe1qpx3eq5a4g/d2+gd18/SJrrXpWl/MNtBrfTHb6B4TeV1n9U/yuvJz6rBnNYvyOrLuCP0596Pdsq4CNfo89/RFnkdu4/nv1d+68+Z17loJJ7xmXf0K3iesD5T3BW2VbTF5/rR/xvmuBk9Xj2j2XXvU7ie8XofGU651unjK+6WLp6yTkud0DNf8XtDFuWx1TzvSB7o1R2dH0rWBpk1zXuxqMOw1jp6WvJcZ43Ku4mfGp4zBtCltmLbmd+f9wpo2HGdNG96DXQ23zp57jaf8vfnbeV0yztEPuvmTfkq/5L9XjevhtCPPg36a15nnkec0e744Nh/tnvkzxs3mxazdOHvOyJpwvOY53zGebv2Z30wbEREREREREZEN4ksbEREREREREZENsil5FFOPMm1xlg6aKaVMk2PKYaY5kWzr9cEHHyzj27dvT4/RtX/cO5mWx9SyWap9VS9jYCocU0C7ltJMmcu0RB4/0w3Zyo/t+vJ8OwkJ6eQ3x5DC2Plipi1+8skny3iW6l012pG+mKms/Ew/r5r7YqY3HrM8KuG9zvs55X3clv5Bn+vkUbQN419KZ0jKZeh/PI8ulTzb5M72S7rU4q3F6Fn86tLA10prcj9+F+2d/sb9GEOrxpjA9PG8t2j/tTY4hhhKutamGf8YG7kWSek3/47+l/MnP2cr7/fff38ZMw28k2RcZPI60Mc4D1aN17Nr/ctj0P/yePy7jz76aNjGz5S55bq5i+3Hxiz25+/meibjH68l412uGznX0mdT2krSNvRNnkfG7lkb8mRPc1/HWjvyOuV6hNezswnnWv5NxujuWYOfeW91tjp0DbNHmRvPs5NH5fMd4yGf03LdyPuAPpXPLbNjV82fM/IY9MXzuv5m2oiIiIiIiIiIbBBf2oiIiIiIiIiIbBBf2oiIiIiIiIiIbJBzrWmT3AvdN7dRo5h6Uh6fOu+qqhs3bizjri3mMWuEU0vJOgWs/ZL1K0jq5Gf1MdI21IbThqnPp8740qVLwzaeI/WQWfeDmsVs58ZzPGZbnwR/b9qR2tPOPvRFbut88datW8O269evL2PqUDtf3IvWdy3pi9Rod+2b+Xd5b/Oa81p2tcQYC7MOGM/j8uXLwzbWp2ENpK5mSldLgjrotPWWa6N0dcI63T1/e/oYmdUCqxrtzeuc8yxr1WRsf/jhh0/c1tUfytgxa5nZafePAd6nWYNmVmOjq9PAGJyxkJ9v3rw5bHv33XeXMX047XRsMfRuzO63vA6zmolVYy2ErrU91zcffvjhMs4YwHP6+OOPp8dYu849tC7fXu6FWS2QPH/GU17HqrEu36yGV0fW5eB90K0vu3pk/Ny1cCd7ip95rl2bb8K4ybVJ1bge4fyU30V70Z/TF3mMrKPKz906lDH2IsXb7l6kT2TN2Fnt2mzlzecCjnONyu/KZ37Ok7yXbPktIiIiIiIiIiKr8KWNiIiIiIiIiMgGOXN51Nq0va41GFPLMr17lkbftRPLFGK2/+paFV+klLZZOny2v+tkNTNZQ+43k05lahpTWa9evTpse/zxx5cx5RppM8pvMt15di9dNKlU1w68S/lkWi//ppPD0R5VYxoj7ZNpwjzHPbe7PIlMEaYMiuPcj9c8r9esdWEn3WAMTR+g1Im+V1X15JNPLuMrV64s44wdlApkaiu/m/dZnkfXXvK86WRuM5smOQetlT3R57gt7xm2Z+e4apSgMkU5U4h5z6xNL96are4nXTydrXOq5lKXjKdsB/3WW28N2yj9pk+tlV0cK/TNtfKMvGac72by06rRP2i7tCO/u1vncpzr3LyHyFo50VbpZDWcW06zRp3JlPIYvM6dRJkxPiU3eczPSWlrJ9Of3bdbkHGspZMN87kjrx/v0/QP2q7zYf4d/TnLKFCew+eOqrEFNW2X8YH+3cXbruX3MXPoMz8/c5x+RBkp15pV45zJtWfG7rWt2tf+zSGYaSMiIiIiIiIiskF8aSMiIiIiIiIiskHORB41SxvKtLUuLXW2LY/N9ESmEjI1qmpME86ONZRrdN2FOvaebprXf9bpJFM8+XdMG83PtE2moDEtjsfL1NNnn312GT///PPDtqeffnoZM7WR0req0fZZUXyWJpf3wV7sm3RpmJ0vzv6uk9TxmlGaWFX13nvvnTiuGjtyHOqLe6S7rrPU7y5mpu/w79Z2HeJ9np2F6IsvvPDCsI2+SFlNdu6gL6ZklRX993ofdPIoSk4zbrJzSaZtZyeFz0n78JiMr9kJ46mnnlrGzzzzzLDtkUceWca8L1LKxlTjnHd5r80keid93htdN6wunnbQZ3n9s7PQG2+8ceI49+WctvfrfVo6+/A6Z9zk5+yaNpOQpz/T5+jb+V1rO9FQqpMSqz3Fx9PSrVl4LfP6U86S8ZSfu66jlANTLvPQQw9Nj9fJIumLGTMpj+r8dK8+3K1vaMe0Ae3IcdVou9w2+y76bNqRcx+7KOYxuDbJdS7turb76V7nxdn6tXuu7J4lZ9LHqrmcPKWiXJekPGpP60szbURERERERERENogvbURERERERERENogvbURERERERERENsiZt/zuamXMWqBWzVu/da3ZqKVLnSg1bdS6VY2thQ9pmVi1H+3hIfCadK3rOr0wazFkXQXqSalLzP2uXbu2jNlWuGrUp1Jb+uqrrw77vfLKK8v47bffHrax5kbXdnFPtj5Ea5r1MWa+mNrxmZ6b+tGqsZZJ54tr2ySmPfZYX4p0sYXa6E4nzdoJVaN/0E5pQ2r3GUOzlgrr1rAuSh6T9U9ee+21Yb/f/e53y7irxdFpjtM3t0RXR4Nka1Pa7sEHHxy2Xb58eRnPakhVjTGb908e79FHH13G2fKbcP58/fXXh23vvvvuMk5fZ82Nta3b98LaGmGpyZ/p+jPuck3Emjapz+c8lnXaeP23qNc/L9bWtOla/9KXOKZPVVVduXJlGdN/83iffvrpMk5bcV3EmMp1StVh8XAvvrfWx3JO49yVtUsY82i3bj9uy3jKc6Q9q8b6ipzfstYbbZproFldsD3VQllbXyprSNGOWWeGtqNN8hiz2lNpb/p9zs98vmAsvn79+rDfrFZq1RjPadMtr2fW0s19szpgVfP6YVmjiHGTf5PXmNc/1yV8zpjVcUzO65nfTBsRERERERERkQ3iSxsRERERERERkQ1y5vIo0rXu6tokMt3t0qVLw35M5+fxsy0pU6VSOkWZAdPT1rbmPAYy1ZQQ+o0AAAn4SURBVGvWkjDbjc5ap1WNaWy0b6bFMf2NKfQp8WAKY0pDmCL+5ptvLuOXX3552O9Xv/rVMmZaf9XFavl9mrRF2oc27XyRpC/yHup8kdf5IvlipsjOfDHb2fO6pm1mUsVM/eV+TNNO6WPni++8884ypo/9+te/Hvajb7711lvDNsbrvfpinhtT3fmbsm0v42geYyadSmkN/47pyjnP0tcznZ9SG9rql7/85bAfbZfzA2UdnTxqj3Ry0076TRt2Kf/cxvs+JaWMA0z7rppLTNdK96q27WP3gm5epH+kXJv+xxbBV69eHfajPIpxM6850/SzfTBtzNiYsYP27qSke/S/vA9nbZPzujI2dvPYY489toxT4kZ5FP0y4y5jaNqQcyHXqOnPPEbOrTP77slHOzvy92Xs4t9lTOW6tJMgMvZyLszj8ZxSWkMp929+85tlnLJhSqdybp3Jo/YkcyMzqXCuL7tnftpqFlurRp+l/2XLb/pfPoMwbq4tw9BxP+1kpo2IiIiIiIiIyAbxpY2IiIiIiIiIyAbxpY2IiIiIiIiIyAY585o2XVvMrpU362NQ00a9YtWofaOuLDWE1LTltpnWvtMX7lV7SHjOqXFmPYvbt28v46w9QV1iHoOaQtqzq59DfWfaifUSqAmtGnXB1Jxmy2/Wvum0xDPN6Unnv0dOo/um7eh/1IBXfbFt4udk3Rp+zjpIM/87zTXfu306X+Q9S1181Vh/gTGtaqyzQE1+6nnpizxG1k6gb+d3MV7Q39IXGUvW+uKeWkXnufJ30Adu3bo17Mc5LetEseYGbZc+y/14HllnYXZvVY1xlO3Zf//73w/7vffee6uOz/Pdkx3X0NW0ydpss3j6xBNPDPuxLhXjdWr36Svpi+Qi1QW7G7zfOL/nXN+1SafNGXvT3rwXGF/TjmwZzLhZVXXjxo1lzPiacbmr/3Uv5tbzpKuF0tUI43XOa871R2drbuPxs+U664DlfPfb3/52GTO2svZJ1ViXqqtps9caRWnH2Toja8lwfuJ1rhrXnqw9lbVq6LP83oyb/G7W6KuqeuWVV5Yx6/TlsxHPN+87rgW6eXEvrI0n3TP/rF5m1rThfrxe+ZzR1bGl33Yt17cQJ820ERERERERERHZIL60ERERERERERHZIGcij5qlFHXpRZlezLQ2pptmuhvTB5k2mulzTEHP9l9MkTzmFO6OLpWfaWaZQsq0M6b3Vo3ymVn77+488pyYwpjphjxH2j5T/pm+n2m0a9vw7Ym159616WMaI1NPM72RvsPrnK2pmQ7ctfxe2wpxz/b5nE6qOItxKW1iXMuUXsowupbfs3PKa0zbpMSNvkhbZzvozhdn98Ge0oe7NHDaivKi3C/9g7IJthJOGRXljjyPTLfvWshTksF2tV1MzXuBv2WvdjwEytMonakabcU08Gz5zWNwvmMMqOr9aLaeOfZ4mnTSmk4GSv/I+55xlPNitoGm7TjPMk5Wjb6d7YMZzxlH83z5W3J+6Oy/B7o5iHEtZTW0R7Z05zEZu9LH1koybt68uYxTvszPjPn5XbRpynaOYV2asZ/3KX0l5yOS5RE4n3LeytbtfJbsymnQxzj3VY1+2s2LvDdSRre2JMce6eaSbt6nb+b7AEL5IH0l7xc+j+bak/Y+RFJ60uf7hZk2IiIiIiIiIiIbxJc2IiIiIiIiIiIb5My7RzG9KFOPmOKWadVMgWL6YKb9M+2JaU4pj2I6YlZrn6Ujdl2DjiGNjaytzJ+pnEwBzDRCStmY+pZdLLiNco2UbvDvug5UHGdaIrdl+vAsff9YbM3flL+d931KKOgftDdTgXMbr3t2x2F6aaYtKlU83Bc7eRRTxJl6mmmo/Nz5In024yTt1vki43+XBr5XX1ybzp92ZOouU72rql5++eVlzPia9uE22j5jL69tppzzPBgDOilb2nGvHU7W0P22tfcs7ZFrIK5heM1znuXn7N61VvK7J7+6F8zmwrQBfTP9g9ea8ZZSmqpRDt7Zm7KetCM/0y/Xyi6qjs/Gs9jFZ4equSy1anwuoMSN46pR4si5L+Md/TS/i+ujmcQjj5nrtGN8Bpmtb1LmRn9JKRI7GlK2n6UYZnNm9zzRSSZnc2Qe49ikislM3p/+wXiV0kI+C9Df8tpxvuM1Trkp1075PmBtp+At2MlMGxERERERERGRDeJLGxERERERERGRDeJLGxERERERERGRDXLmNW2oA0ttGvVoqW8j1Jylxp8aRe6XLdxmetKqUWd3UetoJLP6A2vrNJy0793+vaqvt7G2JVxXa6BrO7sF/eL9hL8ptZu0XfoHrzX/LnXFrJ3B43X6/PR7+uJF0ud33AtfnP1dHoM1Fzg+jS/y7/TF/zHzv86OWbOCvrP2utBW2fKWtuvq3fCc1tYCO+nzsdLF01yLsE4Ka2pkjRPq+hmTs0YY43AXT4+hXfD9oKvB0K1faVfWTOjqhHG8tg151bz+17HXyuhYW0ejq9vFujPdmpKxce06tKv11tmwi6fHyCF2zPUN6990tS9nvpjzYs6FhPY51I7H7JuHPvPzenG+y/qMrEXEv8n5k+3jOa6aP2d0dWzPCzNtREREREREREQ2iC9tREREREREREQ2yJnLo0imGjEtqdvGNCpKMKrmqfhdemmmQM3ShreQGrUF1l6HQ/eb2bBLH157/K616UVue9r5W6YS0nfW+uIshTSPd9FTg09Ld49yW6b3flkfPjS9V188me63H9IuuoupnSxy7Tl2MjfnzD6eptyU8Y9SqWxPy/R9Hi8lc7NU76rjbrl+P8hrxPkpbcxtncxwJrXo/L6TGXbxQf/r7/POhjNp8Emf1xx/rX0vknTmbhxqR9LZqrPx2u+a+d9FlQZ3HPrMz2eQTrrWSbE62dOenvnNtBERERERERER2SC+tBERERERERER2SC+tBERERERERER2SDnWtMmoX4stdj5Wc6PQzV/M83o2nobp/mue3GMi8xaX/zss8/O7Jzki6z1xbV67W6/Q/1eXzwdh2rf12r317K2XpLcnbXxNNuBy7a4FzWBDq2dIYeR15H1K9bW9Dokft7tPOR0HGrHQ2x3mhqA2vVw1s6L2b77ImOmjYiIiIiIiIjIBvGljYiIiIiIiIjIBjmtPOrjqnrnfpyItFy7h8c6NxsqkzgOO15wdmnDL9vi+wjZpR07LmAK99HZ8IJyoeyoL96VzdnwSG12Ekdnx0Nst3N7H50NLygn2vGBnd+cIiIiIiIiIiJHifIoEREREREREZEN4ksbEREREREREZEN4ksbEREREREREZEN4ksbEREREREREZEN4ksbEREREREREZEN4ksbEREREREREZEN4ksbEREREREREZEN4ksbEREREREREZEN4ksbEREREREREZEN8n+ZGg0wefAhdAAAAABJRU5ErkJggg==\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: 4.92 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(1*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.0508, PSNR (dB) = 12.94\n",
      "MSE = 0.0872, PSNR (dB) = 10.59\n",
      "MSE = 0.0311, PSNR (dB) = 15.07\n",
      "MSE = 0.0618, PSNR (dB) = 12.09\n",
      "MSE = 0.0562, PSNR (dB) = 12.50\n",
      "MSE = 0.0285, PSNR (dB) = 15.45\n",
      "MSE = 0.0581, PSNR (dB) = 12.36\n",
      "MSE = 0.0592, PSNR (dB) = 12.28\n",
      "MSE = 0.0792, PSNR (dB) = 11.01\n",
      "MSE = 0.0644, PSNR (dB) = 11.91\n",
      "Wall time: 8.01 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": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1978\n",
      "0.0502\n",
      "0.1803\n",
      "0.2976\n",
      "0.1545\n",
      "0.2693\n",
      "0.1109\n",
      "0.1325\n",
      "0.1272\n",
      "0.2645\n",
      "Wall time: 273 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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3461\n",
      "nan\n",
      "0.7290\n",
      "0.4390\n",
      "0.1995\n",
      "0.7951\n",
      "0.4179\n",
      "0.3019\n",
      "0.1336\n",
      "0.5610\n",
      "Wall time: 170 ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\public\\fyp-code\\plaidml\\.venv\\lib\\site-packages\\ipykernel_launcher.py:145: RuntimeWarning: invalid value encountered in power\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/tensorflow/models/blob/master/research/compression/image_encoder/msssim.py\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from scipy.ndimage.filters import convolve\n",
    "\n",
    "def _FSpecialGauss(size, sigma):\n",
    "  \"\"\"Function to mimic the 'fspecial' gaussian MATLAB function.\"\"\"\n",
    "  radius = size // 2\n",
    "  offset = 0.0\n",
    "  start, stop = -radius, radius + 1\n",
    "  if size % 2 == 0:\n",
    "    offset = 0.5\n",
    "    stop -= 1\n",
    "  x, y = np.mgrid[offset + start:stop, offset + start:stop]\n",
    "  assert len(x) == size\n",
    "  g = np.exp(-((x**2 + y**2)/(2.0 * sigma**2)))\n",
    "  return g / g.sum()\n",
    "\n",
    "\n",
    "def _SSIMForMultiScale(img1, img2, max_val=255, filter_size=11,\n",
    "                       filter_sigma=1.5, k1=0.01, k2=0.03):\n",
    "  \"\"\"Return the Structural Similarity Map between `img1` and `img2`.\n",
    "  This function attempts to match the functionality of ssim_index_new.m by\n",
    "  Zhou Wang: http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "  Returns:\n",
    "    Pair containing the mean SSIM and contrast sensitivity between `img1` and\n",
    "    `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  img1 = img1.astype(np.float64)\n",
    "  img2 = img2.astype(np.float64)\n",
    "  _, height, width, _ = img1.shape\n",
    "\n",
    "  # Filter size can't be larger than height or width of images.\n",
    "  size = min(filter_size, height, width)\n",
    "\n",
    "  # Scale down sigma if a smaller filter size is used.\n",
    "  sigma = size * filter_sigma / filter_size if filter_size else 0\n",
    "\n",
    "  if filter_size:\n",
    "    window = np.reshape(_FSpecialGauss(size, sigma), (1, size, size, 1))\n",
    "    mu1 = signal.fftconvolve(img1, window, mode='valid')\n",
    "    mu2 = signal.fftconvolve(img2, window, mode='valid')\n",
    "    sigma11 = signal.fftconvolve(img1 * img1, window, mode='valid')\n",
    "    sigma22 = signal.fftconvolve(img2 * img2, window, mode='valid')\n",
    "    sigma12 = signal.fftconvolve(img1 * img2, window, mode='valid')\n",
    "  else:\n",
    "    # Empty blur kernel so no need to convolve.\n",
    "    mu1, mu2 = img1, img2\n",
    "    sigma11 = img1 * img1\n",
    "    sigma22 = img2 * img2\n",
    "    sigma12 = img1 * img2\n",
    "\n",
    "  mu11 = mu1 * mu1\n",
    "  mu22 = mu2 * mu2\n",
    "  mu12 = mu1 * mu2\n",
    "  sigma11 -= mu11\n",
    "  sigma22 -= mu22\n",
    "  sigma12 -= mu12\n",
    "\n",
    "  # Calculate intermediate values used by both ssim and cs_map.\n",
    "  c1 = (k1 * max_val) ** 2\n",
    "  c2 = (k2 * max_val) ** 2\n",
    "  v1 = 2.0 * sigma12 + c2\n",
    "  v2 = sigma11 + sigma22 + c2\n",
    "  ssim = np.mean((((2.0 * mu12 + c1) * v1) / ((mu11 + mu22 + c1) * v2)))\n",
    "  cs = np.mean(v1 / v2)\n",
    "  return ssim, cs\n",
    "\n",
    "\n",
    "def MultiScaleSSIM(img1, img2, max_val=255, filter_size=11, filter_sigma=1.5,\n",
    "                   k1=0.01, k2=0.03, weights=None):\n",
    "  \"\"\"Return the MS-SSIM score between `img1` and `img2`.\n",
    "  This function implements Multi-Scale Structural Similarity (MS-SSIM) Image\n",
    "  Quality Assessment according to Zhou Wang's paper, \"Multi-scale structural\n",
    "  similarity for image quality assessment\" (2003).\n",
    "  Link: https://ece.uwaterloo.ca/~z70wang/publications/msssim.pdf\n",
    "  Author's MATLAB implementation:\n",
    "  http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "    weights: List of weights for each level; if none, use five levels and the\n",
    "      weights from the original paper.\n",
    "  Returns:\n",
    "    MS-SSIM score between `img1` and `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  # Note: default weights don't sum to 1.0 but do match the paper / matlab code.\n",
    "  weights = np.array(weights if weights else\n",
    "                     [0.0448, 0.2856, 0.3001, 0.2363, 0.1333])\n",
    "  levels = weights.size\n",
    "  downsample_filter = np.ones((1, 2, 2, 1)) / 4.0\n",
    "  im1, im2 = [x.astype(np.float64) for x in [img1, img2]]\n",
    "  mssim = np.array([])\n",
    "  mcs = np.array([])\n",
    "  for _ in range(levels):\n",
    "    ssim, cs = _SSIMForMultiScale(\n",
    "        im1, im2, max_val=max_val, filter_size=filter_size,\n",
    "        filter_sigma=filter_sigma, k1=k1, k2=k2)\n",
    "    mssim = np.append(mssim, ssim)\n",
    "    mcs = np.append(mcs, cs)\n",
    "    filtered = [convolve(im, downsample_filter, mode='reflect')\n",
    "                for im in [im1, im2]]\n",
    "    im1, im2 = [x[:, ::2, ::2, :] for x in filtered]\n",
    "  return (np.prod(mcs[0:levels-1] ** weights[0:levels-1]) *\n",
    "          (mssim[levels-1] ** weights[levels-1]))\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = MultiScaleSSIM(original, contrast, max_val=1.0)\n",
    "    print(f'{a:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 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
}