Simple Autoencoders in Keras - 784-2-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 = 2"
   ]
  },
  {
   "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": [
      "1595258335\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 33:29 - loss: 0.6928WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (9.064999). Check your callbacks.\n",
      "WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.892000). Check your callbacks.\n",
      "WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.446502). Check your callbacks.\n",
      "235/235 [==============================] - 61s 258ms/step - loss: 0.3873 - val_loss: 0.2807\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 29s 121ms/step - loss: 0.2748 - val_loss: 0.2684\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 32s 134ms/step - loss: 0.2618 - val_loss: 0.2537\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 29s 122ms/step - loss: 0.2512 - val_loss: 0.2484\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 27s 114ms/step - loss: 0.2476 - val_loss: 0.2458\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 28s 121ms/step - loss: 0.2454 - val_loss: 0.2437\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 17s 70ms/step - loss: 0.2281 - val_loss: 0.2272\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 16s 66ms/step - loss: 0.2280 - val_loss: 0.2271\n",
      "1595258335\n",
      "Wall time: 20min 34s\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": "iVBORw0KGgoAAAANSUhEUgAABG0AAADPCAYAAABP9SwxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAGlUlEQVR4nO3b0WncQBRA0Z2QEuzvbBHqvwK5B+c77mFSQLzghXF0Jc75FoPg8VbishpzzhsAAAAALT+OvgEAAAAA/iXaAAAAAASJNgAAAABBog0AAABAkGgDAAAAECTaAAAAAAT9fObil5eXeb/fv+lWeOT9/f328fExVpxlhsd5e3v7mHO+rjhrjDFXnMPz5px28eRW7qI5HmPlc9Hv6aHs4sl5R70Gz8Xzs4vX8GgXn4o29/v9tu/7urviS7ZtW3aWGR5njPH76Hugwy4eZ+UumuMxVj4XOZRdPDnvqNfguXh+dvEaHu2iz6MAAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIEm0AAAAAgkQbAAAAgCDRBgAAACBItAEAAAAIGnPOr188xtcvZqk551hxjhke6m3Oua04aNu2ue/7iqN4wrZtt33f7eL5LdtFczzOquei39PjjDHs4gV4R70Eu3gBdvESPt1F/7QBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAIJEGwAAAIAg0QYAAAAgSLQBAAAACBJtAAAAAILGnPPrF4/x53a7/f6+2+GBX3PO1xUHmeGhzPH8zPAazPH8zPAazPH8zPAazPH8zPAaPp3jU9EGAAAAgP/D51EAAAAAQaINAAAAQJBoAwAAABAk2gAAAAAEiTYAAAAAQaINAAAAQJBoAwAAABAk2gAAAAAEiTYAAAAAQX8BgZfUlIzx7HsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dW7MdVdWGByqeQQ4SjoEQAoFIIHKQG2/EK6usUv+qV5ZVVlmCUMgpSAIESEhCQjgZxfOZ7+KrtM982XPS2exk91r7ea7m2t27V685eow516rxjnHFxx9/XCIiIiIiIiIisiw+t903ICIiIiIiIiIin8QfbUREREREREREFog/2oiIiIiIiIiILBB/tBERERERERERWSD+aCMiIiIiIiIiskD80UZEREREREREZIF84WJOvuKKK+wPvk18/PHHV2zFdbThtvLhxx9/fMNWXEg7bh/64lqgL64B+uJaoC+uAfriWqAvrgH64lqwoS9e1I82IpebK674X+z5+OOVjx+ntvsGRKSq9EWRpaAviiwDfVFkGWzoi/5oI5+Jz32uVdh9/vOfn3XsC1/436PHH2aqqv773//Oeu/Rjzg8xuuN/uc///lP9zX/bw1+PBIREREREZEVwJo2IiIiIiIiIiILxB9tREREREREREQWiD/aiIiIiIiIiIgsEGvayCxYn4Y1Xb70pS81533xi1/ccFzV1rHhOGENmn/9618bvm9VW3Mm6+Dw2Nxixqy5k//H6//73/9uzuPrufV4djJZw2irsebQpYF2y3pVc206qg1l3ahLB+0z145p05590lZza4jJpWEr4qt2235Gvkjm1vaTzZPz37PN3POq+rYZrYu5v9S+W8dm7Ujc38ilxkwbEREREREREZEF4o82IiIiIiIiIiILRHnUBmwm1X8dmJse+OUvf3kaf+UrX2nOu+qqq6bx17/+9ebYddddN+s88s9//nMaUypVVfWXv/xlGv/pT39qjv35z3+exn/72982HFdV/eMf/5jG2fKbqaj8/CntGqWv7iRGfkMZxqVIPd2KNvHyyfmn3SgfTB+48sorN/yfhD6W/kaZYUoQ81z5JCPZE1/TdrRbVStpzWvQd2iPtBXjdMZs/p++eHH0YuPFSBV7UuGLSeXXblvHXJ8d2Xhkq5GEXDt+kp49UjrP14yhuS6OygD04mnGTL7OdbAXT7Xtxlzq/c3cMgp5TEnxcphbTmO7MNNGRERERERERGSB+KONiIiIiIiIiMgC8UcbEREREREREZEFsmNq2ozqaMxte3oxuu9V1CjyPnNO2Nqb9WhuvPHG5rw777xzGu/fv785dt99903j3bt3T+NrrrmmOY96Ump2WcOmqup3v/vdND5z5kxz7NSpU9P4+PHj0/jtt99uznvvvfem8R//+MfmWK+V90gXmzUiUp+8DvDz97TdVW3tI9Yt+trXvtacR3vnfLHm0N///vdpPKpNNNILz/XndacX41K7zxonrF+Vdaj4mnZPaJu//vWvzTH6d9qXr/mMpK3XndE6Rj+i3apa+zDefvOb32zOY2wfxTLa7ve//31z3kcffTSNWVusar4dd5IvklEdk54vfvWrX23Oow/nPHLOGTNZOy5fj2q92W7600mf7a2ZaUfaOGts8Bq0D9fIfJ1ra8/G616Xb7T3pz04/4yLVW19Ru6B+feqdi1MH2Ns5F72/PnzzXmMr7lm9moyjnx2J9CzcfpR73tNfifhMf5PVRvn6G/5fYJ1N9OOPT9NOxpTN8+oduao5maPUY2wS+lvZtqIiIiIiIiIiCwQf7QREREREREREVkgay2PGrXsY6pxpoEzhS7/j4zanvbau41aoG43I5kP0+tvvvnmaXzw4MHmvEceeWQaHzp0qDlG6dS11147jTPdkDDNLOeO6aUffPBBc4ySKKaoZvrwH/7wh2k8anE7SgkfyUtWlVGaPp8N2i5lT7QxU4hTkkH7MN23qk0Vfvfdd6fxhx9+2JxH+4zSFjme2158HRil6DMWfuMb32jO27Vr1zS+7bbbNvx7VZtOnP7M9x6lD9OmlC1Wtf7NZyIlk4wR62LD0To2V7a6b9++aXzvvfdO49tvv705L9P7Ceeavnjy5MnmPL4+e/Zsc4wyAKaLp99TSrAudrzAqLVzL7ZWtfa9/vrrp/ENN9zQnEdZR66Z9DnaIiVuW93ye91suBE9SUbupbhO0nbps1wnU3LKNY6SQ+5nqlq75jHunxiXU8bD/euq2HEkgeL+PueV69hNN900jbl3raq6++67Nxzzf/L6uVadO3duGh87dmzDcVUbT99///3mGO2RduuxLlKpkdylZ+OUPd1yyy3TmHa84447mvMYb1NixRjLmPrOO+805/E1bZ//1/PLqvVeF7eCuX6fMZl7YI7T1rx+fnfvSf/Thp9V0m+mjYiIiIiIiIjIAvFHGxERERERERGRBbLy8qiR7Imp/inJYJpcpkoRpkCNqr9nNXCmrHKcaXxMC9/utEWmk2XaaK8TECUwVW16fXZDYFoYU3UzXaw3J6N0tOxSkun2c9hsuuGosxSfraV1khrJupgWmGn6TL9nqjf/XtX639VXXz2NM52fvpjPAv2ZNs50/s0wt0p81Wqkoo4q4mc3Idpmz5490/jBBx9szqOUhqnEaeueHDTvi76SPkpJFOWNVVWvvvrqNGYcHslNlyQ9/TRGvsgYkl1MmI7P7nwPP/xwc979998/jSmVSjkU7yPXO8qZaCvKZas+2VmM9KTCaauR9HgVmCuBqmp9k+ts+hhjKPcz9Muqdk3OPQXT8ul/KVUcyaN6rFs83YhR+j3tSD9N/9i7d+80HkkVR37EPSXloim7GEkVydxOqEuOqb0YmvsXzmvuRSiLoVyGNqtqu5/Sbvk9g89EdkTkfTAWsvteVSsbTokbPyefx4vpdLsqvpj+xs8+krlR2pb7G5Z3oC/mczHqcMrvfqOOttzD5BpAH+N6t27r4lYwirvc1+Z3U66fuR9mjOitx1XtnKc/0zcpY8zvKpRJbuY7v5k2IiIiIiIiIiILxB9tREREREREREQWiD/aiIiIiIiIiIgskEXVtBm1cKMGkFrQbLFHXf93vvOdaZyaVOrFU7tPDRr1wqlNO3HixDROvTBboo40wdRHXu6aNqPaCXkv1O9xHrL2BPWGqa+m3agNZK2EqrZ2Cecn75d6w2zNxmucPn16Gmdr8N57VfW1viMN8Kgd+HYwqnMyqlvDuR3VN8p5J702onlP1Jrm9Vh3gc9TXoN+NWr5PYLXXLLOu6fdT500/S3byTJOPvroo9P4gQceaM5jPQZqfdNnGSezVg2fOcbdfOZ4j3wmqlrdOH04a3Fst79dDL2aQ1kLjHVnWEuhqtXkP/TQQ9P4wIEDzXmsu0Ctd85XtqcktB1rN6S2ntpuPhf5mrZLffhoL7Ak3+y1Y09f7Gnm8xj/L9c7wjiZvnLrrbdueF5ek+s4axRVbU0NjFWJp0lvjaxqbZU1Zxi/WA+FtTLy9f79+6cx2wpXtbbKvQnrnjAeZi2OUY2qXlvaPI/zsaSaNnPbPLPORVXVrl27pjF9paqtacO5zOeAc8eaM7n3yFhO6BO8/mhPldA2o7VvVeJpVXt/jIe5X2Dcy7pR9DGui1nT5q677prGfC6ybgr3NOkfXLtGNR7pw9n+nbF47rq4boy+89OHabOqqscee2waP/LII9M460vRN3OfkzVpN7qHqvb7Yn6/feWVV6bxqPYNX1vTRkRERERERERkTfBHGxERERERERGRBXLZ5VFMO8sUNKYiZVohU0+Zwv/973+/Oe/xxx+fxkx3y7ROpvdnGz0eY9pUnsd0rjzWa7+3hFTEC/edaZj8PJmO2/us2daOqYN5jDAFkO3R8hjnkW38qtpnIlvXcp6ZPpx22kzL9blpqHkfl4tRW9Je6vcojT4/E23M+UuZBP2ZqaI5zzyWKedMJRylGvOaS/CxrWQkcWPKMFNzq1rpaKboUzrKMdtBV7X2YBo4paH5OmMHU1vZnjjTV3lePkv0dT4vq5QunL5In+NnSrkL7ZipwZRhsLVpyuF4fa5vtGlVG3tzbinr4bFRa+r0WcoWRjGGLMmf58bTTKvm64y1/HyMp/m5uYehj+X8069SUkC/omw4P9c6x9OqcUzlc55zyz1ISjLuueeeacw9KiWMVVX333//NKYkKlP2mYqf+1f6FWU8uba+88470zhlznwOVyWO9qTBVe3noQ1H0u+UKnKeuVfMNtx8b8pN+Z2jqrVv7m1o71GbZx4bSb9HEqje/yyBuTK30bpImWFV1aFDh6bxt7/97Wmc+yBeg3ObclG28k4b9PbR+X2Fz0J+FvpzxuJVZ/R9pFf2pKrd6/zgBz+Yxj/84Q+b8xh3eb2UJXGvk6VO+J1mVBqCsraMK4y1/JybLdfQY72eDhERERERERGRNcEfbUREREREREREFog/2oiIiIiIiIiILJDLUtNmru6b+uHUhlKX+KMf/WgaZ00b1kWgbo0tuPN1tq/t6exSn0+tXrZwox6Z2rpsF7edLRRTXziqvUNdLbW4qfUlqRvk//FYzj+hPfN+R7VQ+F6c46yVwWMjreFmNcJL04r32rDmc8jnNOeM/8faCjkPWWPlAmkr6opTQ0o/ooY0daKjGgxL03BfLPkMUU9PTf6o9S9bfFe1rS/ZHjrbztJPX3311Wn8q1/9qjnvrbfemsYZ1/fu3TuNafusCTGq68LXvH7adjMtFC8XIzvyM6XfsGZCtvTt1fXKWMm6Xqxl8vbbbzfnsaZN+intNao/NKo70WtRu2SfHdWK6MX3fA4ZJzPW8rOO4livHmDaiTVtRvGUx3ZSPK2aH1Nzz8c4xHWrqmrPnj3TmLUzMvbSv0+dOjWNWROhqt1bpT9fddVV05h2HO3pstYY13g+k6NnYUmMfHHks/w8rCFV1a53HOd5fC44r6OaHaO6UYyh+R1hbuwY7d+X7LM5Lz1fHMW5bN1OX+Q46+ixZukbb7yx4biqrWmTMbVXp29U/yyP9daRVdrfkF7b9qwhxXia9RR//OMfT+Of/OQn0/j2229vzuNccv9y8uTJ5rwjR45M46xtyn0Va+mM2rann54/f37D62dtnc9qQzNtREREREREREQWiD/aiIiIiIiIiIgskMve8nuUwse0OKZ/VrUpUUx3Y7pvVdXRo0en8SuvvDKNmYZa1aY+Zjo6U1vZUjVTE5lynqmtTI/qSXWqtjdt8WLaUzOlayQ3oj3yepSQMS0x6bUIpi2q2jZ/2X6N9mYaW6a0zW1tOrfNd7ZzzZTk7aZnx5EcbnQN/l+mPtImbH/I1sRVrV2z7enZs2en8Sg1cSvS+UdyjSXRe94yfZi+w7aIVe2cU36TUsXDhw9P4yeffHIaP/30081577///jROiRWvz2durtQk4TqxSmngc9P5M4YwHTvXxbT5BZjOXdX6EdOGM4WY8XskX+M9ZaxgbM/U4FWXZIyet1E8Hcmge89sxlO+5p5lJA3IeMpU8s2ui5vxsSX75cXQa0dc1W91n7748ssvT2NKFbPNMP0oZR2UnFIimfbuyX2qWl/n/+WeZUm+OFem3pMeVbWflXKKPMbrjWzN9ShjN302paKcZ75vfqcZ2aYXQ0c+u3Rf7LWIzrnlejRqr819KL+zVVUdP358GlP+/frrrzfnUaqY62KvLXk+q/TntDGPjdq/L912F+j5acrCuH/Jtu2PPvroNGa5lPTnc+fOTWPK9n/+858357355pvTOJ+X733ve9OYe+O8X36HPXHiRHOMzw/3wxmTbfktIiIiIiIiIrKG+KONiIiIiIiIiMgCuSzyqF7aYqZ+8Vh2IGFaOKs0UwJVVfXMM89MY6YvZYoq0xYPHjzYHGO6G1PwUgLFjikffvhhc4xp4dvZIWpEpmnNlQBxnJ+NnzvT2JiCzVRHyqGqWvkMbfPII48051H+wTTjqqpjx45NY6Y2Znop52CufGGULnwxafGXg1Eq3kieMBfaMTttsJsN07mZ9ljV2jtTDinXYWprphzOlbmRVUk1Hd0n5z87Lu3evXsap7SQ6ab0y+yawFR+pn+mv/EamcpP/2YcT/kHn8dRGvjouV2yTefKTDIOMb075VCUyXBeuEZWtb7D1N3sqsj3Tskp34up/tk5kT6b8oNeJ8VRh7rtZq5suLfPSUbrLPc9GU/pR/RtdoLLYyl/o224Zxl1uFi3eFo17srCZ3H0XOYelf7C+aS/VbX7Ru4pU7rB+JgxlftS+mL6G22ce2BKNNgdacmSjLkdkmjPXEu4d0jJEu3N9TTnn6/ZeYZ70qq2w1i+FyUxtNtIsjVXHrUu6+IojvJ7WspdaDuuM7necb9Jv8y4SRhfq9p1kT6b6yJLZuT6zLhMe2f8WZJUcUTPhnn/tGGWTeAelfE0954//elPp/HPfvazaZzfJegr2dGPEn4+O3m/jNfPPvtsc4zPDH17q21mpo2IiIiIiIiIyALxRxsRERERERERkQXijzYiIiIiIiIiIgtkW2vajLRe2WqLukS2Y2Ptkqqq1157bRqzlknqyKlJveuuu5pjbKHJe2fLsKq2XVzqF7e7lskccv6zpR6ZW9NmVLOIml5qUHP+2XKddWyy9hDvI/WLtD11oflcpc64d33+30gPvwS7j/TLoxo+ZHSMc0FtfbZ6ZmvEBx54YBpn+2leI2vVUEPKlqWjluSbZam675Hmm89vzj/96rbbbmuOUXvNlolsDV3VasA5/xk7WKuG9Yuq2no6vI+saUPbs8ZCVasR5ngJ/jaXUVtr2jTjMGubZN0iHuM10j967aizRg7rM2QdJPozbUetflVru6yVQhtzTV+VGgxbEVsT+jA1/lnrjbY5cODANN63b19zHmsspB/R11lDZck1hS4Fc9fwfC5p19w78HWvdlBVG0fpD1nDiPvQ3COxzTc/S7b1po1HtVJo/yXXzRjVO+nZcG5txao2HnKOc/1kDYw9e/ZM42xbzO8ZWcfkzJkz03hUx2RUx6f3fWo0N0tjVAdytC7SVhkruafkGpQ+wLpOfC6yxTtjb9r49ttvn8bcB+W6yPfK++D3Rz4Lq7S/IT0/TVvThlk3in7AejG//OUvm/OeeOKJacz9K+1Z1doma9rwOwltz/Wyqq1jc/To0eYY7Xsp/c1MGxERERERERGRBeKPNiIiIiIiIiIiC+SyyKPI3FS/TAvrtUTLNDOmbTOdO6UD999//zSmdKOqTdOi7OaFF15ozmPa6yqlsV2Y90zDZPphyohIrz1mHssUcaYsMt2Qrd2qWukUU9oyHZnpaO+991732KhlLu8x2ykSvnemvi05tXizae6jVH/OBVO6M0WV7S6Zfpg24Pxl21PKc5jmuu7p+3Phs810+qo2bZvp3FXts8607UzpZaxlTMgWm4y1KWNkOjH/L31l1AKakihKPlbpOZjb8nvUajvlUVzvaNOcW/osn4X0WT4zKY9izKYvZtr6KB26Z+NVsiPZzH1nPO3JTdPHGE+Z3p025F4k2zxznbyUbUlXjdG+tMdIWk3/yH0F/ZstotOOXDOzHS7jwKlTp6ZxSrEYz1Mqt+r+N/f7Q+7NR9Jq7jevu+66aZy+yBhK6RTjZ1W7P6JkP++fMSDjKY9l7NjMc7s0LsW6yDnr7R2q+iUb8nqUfFOaWtXafyRz432kdKcniVpVm5Lec17V+seVV17ZHGPsopSQ46rW7+m/KTflPvTxxx9vjnHvzO+Ohw8fbs578sknp3HKHS/XGmqmjYiIiIiIiIjIAvFHGxERERERERGRBXLZ5VEjmE6W3ZiY9snURKbIVbXpiUxvzOrvTHHbvXt3c4zpaZRHZfcopsCuYnpx3jM/d6aq9ToqjNL38hq0FdPAM42Q88p0xkwvZWXv7B7F/2PacT4vfM6YFlfVl7xluvNISrTk9MbRvc3tHsWU1UxfZYopfTbfl/P+9ttvN8coixxJ73YS/Ox8nlMexe4X2amJvkQfyPRh+jD9KK/H9OHsZkO5KVO/mS5cNe4etW6ymqp+2nDOLbsZZFeLnlQ1YxT/j3LUTCFmJ4yUrdKfGRsznZ/PTKZD92S3q2zHjbiYz8M54dxlPOV+hlK1hH7F7ntVrd/Thjs5niYj24069DB+0Y7pH7QdJRmUv1W1kqi0N9+ba2R2j1qHTjQ9RrKaUZcpzkNK3LgvzVhLGKPplyn95nuP9pd8r1GHwNwDj8oYrAOMS/l9gvEx10zuF7gW5jPD/ROvn+viqBMnfZjy0/xeMyqjsA4ytznk88o5yVInXLsolUq/5J6FvpPnffe7353Ghw4dao7xOXj11VenMeVQVVVvvfXWNE6J2+VivT1eRERERERERGRF8UcbEREREREREZEF4o82IiIiIiIiIiILZFtr2ozagWbrWbb+JWzxVdVq06gRzlZ81ChmO1zq1l577bVpnO0z11kjvNnWgtQIj1r0UYM6qgfE2hbZDvrFF1+cxsePH2+OUSs5qoHE66cenMc4Tm0y7z91jpzHVdKqjp6FuTrqnOsLsGZRVVubKGvacK45z3kPO6kmAz874x9rx1S1+vqMVaxjM2qlTn/m9aj9r2rr52TrWl6TfpS+Qu15avf5mv+3irXELsBnlvM8qhGTjOoAEcZiPj+jFuL5LPA54ftmXYjR/c6th7bOZKzqxa78O+Mpj6XduU4ytla19U9G+5dVXbe2mowvjD1Zd5E1GUYtnGlH7k2yhhTr2GStFNqVdTRo36pxTZtVr6MxWvNHNW14LG3D16P54frHeJd1TPi9JetLpa0ukPVURm2RR+3Ayar6Mz9frjM8lnsJ7hvpA/nZWfeENs09zKieH69Pm45ay2fNud66uEq2mkPGIH7neuONN5pjrM3G+UmfpW1ow1tuuaU577HHHpvGWT+Mz8uxY8emMevbVPVrbF5OzLQREREREREREVkg/mgjIiIiIiIiIrJAtlUelelFTBlLedSpU6emMdOB2W6tqpUIMPWUbfmq2tTTTCs8c+bMNH799dencaYzrps8iumGmZbHORq14WPqYNrm+uuvn8ZMaUvb8P94HymBeu6556Zxyudoe44zZY6fOZ85pj8zjS9Tzket6teR3rOQaaNMA++1KK1q7Urfq2pTzulv65Y2ejEwPZTpoDn/nK+UG50/f37DcZ7HtFGmIGeqMl8zXbiq9QnKdDJ+Mq6nP589e3bDe1pleRRjz0hyyGc94wufBc572qAnLcy264yBaR+m6c9tX5r3sc4tiDdLT9aQaeCcf85dPhOMofSbqta+o1a4Ozm+kpQ48PlNOTVjFmNlyl3oL5S3ZgxgTE0/pV0pj2Isr2oljSNJxjowVyLN80ZS1Ln7XNo61y3KK958883mWE+2k+s418xcd9eR3rqY6wXnPX2xd16uR/Rvvlf6LP8vy2TQdtzbZhkAXiNLLOyUdTFlbPSXZ599tjnGUgm0R8Ytzhe/c2RrdrYGTx/junj48OFpnN9HuEfdrjXSTBsRERERERERkQXijzYiIiIiIiIiIgvEH21ERERERERERBbIolp+U6uWGl7qdnksW5ZSU8g2tNkyka3eWDejqurll1+exqdPn57GqYdcN933qK5Cr4Vltv6lPW6++ebmGDWGrCmUrYpZ74b2pF2q2tbsqdfutX679957m/P4XKR+ndpV6pGzxeq617EZtfymNjTbJPZaTtOXq9q6UXmM2l/aeFRHY92hvp7j1FCzvkHqpE+ePDmN6UfZcp2aXsZJ1mLY6Pqk15Y0/e3DDz+cxmy7WNXG4VGL8iUz8iPWKsjPxHjDOapqYw/nM+eWWnJef6TdT1IH3nsv1kXKGhu9+10lO241vRob1OdX9WvaZI2FEydOTONsM0zbjOLpTqNXDyXjGv0jfZG2474l96hc0xi/b7311uY8+n3uh1lX7/333+/eE+8j90jr5nO04agVNu2UtWp4LmMyvy9UtXNH/8v185VXXpnGH3zwQXOM1+R7ZZ2dUbvuXqvoVWbufoF+NKrdRR9Of6ZPcN4zHjL25jrIe+Q6m3VrWL8l/ZR7mnWOxTn/3FNmy2+uXazrlP5Be+zdu3ca8ztmVb/GZlW7Bz5y5Mg0zjqnS6g3ZKaNiIiIiIiIiMgC8UcbEREREREREZEFsq3yqITpfZmG1Esfy5RPphwyHSpTVPleTI2qqnr++eenMdPYlpAadbnIlFKmbzJVLWVnu3btmsaZ7nvnnXdueCxtQ1tTxsFxVWv7lGuwvRslUQ888EBzXrYbJ5SGMO0xW34zjW9dUht7bb2r2meBMrS0I32RqaIpmaBdU6rIeR+1GV5ncv6ZDspU0UzH5XOZqcV8tvk8pzytl2bMGFDV+tG1117bPUa7ZYrw0aNHpzHbLla1z8yqympGafq9cVXrA5Q7VLUSCs7tqFVxT0ZQ1coFUu54/fXXT2PG6JRUUSKQz1Ovze0q2fGzks8BX9O3U7rGY5z/jKfcz6SPcf6VR/2P3nqX/sEYmJIlth3mednmluT+tfdeaWPK3CjPSXvzHtfNxnPlpiln4b4k5VF8zb1NyqNoG875SJ6W1+CelbbJtt4j2WsvbubcLDm+jvaXtEfOC2Ngtvzmcz/6Xsl578lPq9pnKPc+vMeRNJjrYspuenuaJdttM+TnYfwbycJH8kH6Eb+PsixGXiP99KWXXprG/D6SpTCWgJk2IiIiIiIiIiILxB9tREREREREREQWyLbKo0ZpwnmMaWyjauBMlbrpppumcaZIsuL+U0891RxjxfdV7VQyh4uZ/7kdLiiNyO5Rd9xxxzSmPCrT3U6dOjWNmUaYqaHXXHPNNN6zZ09z7ODBg9P4vvvum8Z33XVXcx4lAJnGzE4nmdpKRrK5JT8zvY4ZeSzTUulLTCHOtFGmPlKCk/IypipmOuJOlUSRURo4yec3pTSEz2xP+ljVdnajn+7bt685j9LH9HvGZKbys2tYVdUzzzwzjbPrTU8GtGT/SuZKTjPW0N9yHeP/8XoZK5k2TDtmh4Xdu3dPY0pM8xrs5sVxVZtenM8gn7tVst3FMoqt6b+0x8jWtCljaEoDmIafnQ0Zk3eS3PvT6K13o04lczRO/soAAA0JSURBVLtt5vrJPRL9L2WltFXakfsiShDnyi7Wgbl71Iyn3LNyD1nV2oBrX57Ha4467nF/OZJpURY+6nCUx2hTfv5V6q452t9QepR2HEmn+HrUmZH7V9oq1z52IqZMuKr1Mfrp8ePHm/O4p8lnZqfsbUc+m7bp7RXyGtxf7t+/fxqnz9I2o70n19Ml+o2ZNiIiIiIiIiIiC8QfbUREREREREREFog/2oiIiIiIiIiILJDLXtNmVDeFWsbUC/c0wtnq+bbbbpvG1KRmK2Hqvn/96183x959991pvJN03yN9Ied/pNenljhtw9a/1BtmLQ7OOd83Nd9s6XbgwIHm2D333DONWWMjW5Tzs7B+UVWrD2ctjqzTsEp1j3q2m9t2sarVZvNYXoMtqGnTURvg0X2NnsHRM7nq5GfrtZ3N55Lzmu2D6UucV/poVaunp/47a0hRA551cRh7qSV++umnm/OoAc8aDqtax2YzjFq8j2IqazVk/YTe9TKmMlZmy28+X2wrzRpwVa1/Z42qnh3XwaZza8Jl/QXGUNomr8H5Z00+jqtafxu1CF7nmPlpzK05lH7EGJh7Cfom/ZJ70qqqvXv3TmPG0awPyH0F9x9VVWfOnJnG3K/mHmada2XM9bHcv3AtzLoXu3btmsasiZm1v8jIZxlD8zsN12faKdc+nreq9RM3S2//V9XO+9VXX90c47zTT3MfxNdcC2n7qn5tv6o2/nJ/8+abbzbnMX7v1O+VaUP66dwaVbkHYp3Su+++exrn2sfv/E888URz7MiRI9N4iW2+iZk2IiIiIiIiIiILxB9tREREREREREQWyLa2/B6lSo1SiJkelamnbEXL89hWuKrqhRdemMbZmo3piOuYcthjlDrN9E2mnWUqGecuU3V5LlNAKfHI/2MKJNuEV7Upq2xVW9WmPfL6mUrONMWzZ882x15++eVp/MYbb3SvwfvNlLztZm5b75EvjuRRTFvMlE/OO2U2KYcbtZzmaz4/eb8jacKq+3DeP+eBEr4PPvigOY/+kendjI1My8/34lyO0sp5XsrfmJb61FNPTePDhw8357HVYvrROsqj6C+0acoFeSyfbaZ0sxVpphDTxyjxSH+jr58/f745xtTv5557bhofO3asOe+jjz6axiN5zjrYcSvkprRBpt4TPhf0+5RF0lfyvXptcneaVGruszeSZKSciRIK7lW4J61q0/kpR0zbU/ZEOVRV1enTp6cxn4Wl7T8uJSMbjo5xnlP+Rpv2xgmvkW2pudehXLyq3R9xzcz1k/83d11c5djKz8R9I8d5XvoO10LufVL+TVkVJVUpo+K80y+r2u8JL7300jRmi++q8ffKdbHdRsyVBo/ato+k+YcOHZrGlLVlSRRKoH7zm980x5be5puYaSMiIiIiIiIiskD80UZEREREREREZIH4o42IiIiIiIiIyAJZVE0b6ttSG8q2bdQLP/jgg815bPVMnWPqENmyNDXhS9e0bRXZDpKfO+uT8Bi1mTl3586dm8bZEpNQW5p1cajppV479ce8/9Qv0r6sQZPaV+qKUzf+2muvTWPWPcr6SEuuaTOqUTL6O31zVNOGutN8Zvic0MZ53qjeA+9jp/hlkn7K542+krW5OJfpY9R8j2qccP6pKc/6Oaxbk77I+3r++eenccZk3uM6tqrNz8R6B4yjObeMqdmim5p82i5jJV+znljaivVosrbC0aNHpzFr2oy0+6N1ZN38eVRXaxRPaTeO8xpcW0ZrDuNprne9OjajdqvrZqeNoG/ymR215k0f436H+xvWuqlqfZYxNWtIsY7eiRMnmmOMndzDrGPc7JHPJW01qqFH38k6M6xtwX3eqEYV3zd9kdfPtu2nTp2axmwPzXhf1cbrdaxpk/dK2/Gz87tAVTuf6WN8Tfvkd05CP8p1kd8haLeqqhdffHEa8ztD3i9tt0r2+ayMaqUyhmYdIdaRYu2vhx56qDnvW9/61jTm2pr7y9/+9rfTOPcsS/veNsJMGxERERERERGRBeKPNiIiIiIiIiIiC2Rb5VGZIsZU3pTWsLU306EOHjzYnHfjjTdOY6ZFZqoUUx9TOkB2Upow0wiz/RrTx0bpoL22pFVtui9T4TKll9dkOmM+E0xZzJRzQrvle/F+KQ2oalMkRy1W+Zytanry6NkeyTqY/puyDs4t0/4zXZnznOnKTGXm/41kF+sO04eZIpx24nOa6aBsg8n23eljTF/l+zKVuKrvK1Vt7OV9MPW5apw+vA72HaWB01Yp0+S8pH9QwkTZcKaLUwLH+0g78nnKNfPkyZPTmPLTjJujeLgOdtwMIykS54TPRNqaexbG0/Qjvs728bw+Y+g6+tvFwM9LfxvFuauuuqo5Rj+g7TIe0nb0sfQ3Si2OHTvWHONay/3rTrJjfjbGGsaglLqMJDL0F9ojW0WzfMNIws1nIvdHZ8+encaURKX8ns/gXHnUKpH3zc/IWJayMf5f7im5nnI++f2wqvVhzmXu8blvoRS8qpUu8h4z9u7UeDv6/sXv/Gy5XlW1e/fuaXzgwIFpzBbfVW2bbz4HuY+iv+WaSZb+nd9MGxERERERERGRBeKPNiIiIiIiIiIiC+Syy6NGqVJ8nV1MmNZ27733TuN9+/Y157HiNNPisjI/01cz5XCUyrzOMHV69Lk5X6PK/Dnnp0+fvuh7YvrcqCNKyrmYssrPMkqNHckX+DnzM4/SHleFuRX88xjTSDNtlJK17GJCmEaa6eg8NpJHrUOa8Fz4WZkan2nVjHGZKso4yVTvlBnSbvSxUUpzpgXTpnyW0o92kg2r2meYdkw5Ra9bWFXVkSNHpjElp9l9kXbNOEro23zfqtbXOU6fXXeZG+l1xhh1tsl57fkz5b9Vrd3ol7l/oT1GclPaelVlvVtFTx6VsYzdhXLeGdt4Xnab6UnDU2ZIX09pDfdWjKPr5l8jRj5G0t/43Gespa16/lbVxlN+V8nzaN98lnhfHOd5o86bq9oxakQvHqatOE+593n99denMbu1cd9T1dqYNs15pn1GHRe5LmbZjZ22v7lAL7ZWtb6Ye5a9e/dO44cffnga79+/vzmP/0d59+h7WtL7vrhEzLQREREREREREVkg/mgjIiIiIiIiIrJA/NFGRERERERERGSBbGvL79RRU9+WesBefZHUAVPDzfZrWdOBmuN8r3WoUbIZ+LlTD9hrqT1qu7gV9Qx4H6lNpg4xa9Ww/kZv/GnX6Gkx+SxWjbWSS6Nng1Eb2vx81BKPat/0ajyMnotRnaveeCfDeRjVtkjfofZ3ri6evjOqfZPPUu8eR360k+Ju1TiG8FjWO2Bti5EdaRPacVSrIZlrRx7bSX7aq2+Tx9LW3H+MamCQkd/zWL4XX2/1Wr0uzF378hjXP+43R3sOHhvFzdyj8jX3SOlvO8mOo30E6bUGr2ptOHfuevasGtfy691T+vMonq67fUfzwj1N+sdm1kXaKtfBUT0/xoHe99Q8b93t1iOfX9owaxbxuz33q++++25zHm01Oo/1hnLPskr1hsy0ERERERERERFZIP5oIyIiIiIiIiKyQLZVHpVpSExxyxaHzz///DRm6ilbgSdMr8qWcGynmS1LmTq19FSpS0Wm/mba5wVGaeCXmlE6bK9t+0gClZ+5l+qYaZqryty00ZFNR2mFvMbcVPy50inT+T/JKHV6K/x05G+jtrM9++4k6czFMLLj3Dmb64sZ8zZjx1Eb2p3E3Dg2khTPTaEf2XAkI+3F0J1qs09j9GyP5P29/cfoWJ7H955rY2Pq/zOah9Gc9/xgs3ubudL5uRLxneyno3Ux9/W980b0ZE6fdr2e7XaalG0OOQeUm77zzjvNsV/84hfT+PDhw9P4hhtuaM7jd1OWR0m51XvvvbfheVXtd7ql28lMGxERERERERGRBeKPNiIiIiIiIiIiC8QfbUREREREREREFsi21rRJqAHM9sF8ffbs2ct2T/L/sNYBW+ONWm0nc7XXvVocqT+e2/Kb/zeqW5O1akZ1G9aNkY46mavJn8vcWhBycWx1zQrr0VweLoUP9HxzFNf0vc2z2XiabWI3w8hu2vSzsdX1Y0ZrprbaGi7GTj17jGrfzLWTfvnZmFtL5mLo1Sb6tPeWrYHzmrVl+frcuXOX7Z6Wjpk2IiIiIiIiIiILxB9tREREREREREQWyMXKoz6sqlOX4kZkyB1beK3PbMN1aXm9DSzKjptlh7eKXQsbys6y45r66VrYcIfH06o1seNmWCN7r4UN57b8XmPWwo5zWVO77igbrjEb2vGKNX1oRURERERERERWGuVRIiIiIiIiIiILxB9tREREREREREQWiD/aiIiIiIiIiIgsEH+0ERERERERERFZIP5oIyIiIiIiIiKyQPzRRkRERERERERkgfijjYiIiIiIiIjIAvFHGxERERERERGRBeKPNiIiIiIiIiIiC+T/AI8HfsiKgwzpAAAAAElFTkSuQmCC\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: 14.1 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(2*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.0388, PSNR (dB) = 14.11\n",
      "MSE = 0.0603, PSNR (dB) = 12.20\n",
      "MSE = 0.0183, PSNR (dB) = 17.37\n",
      "MSE = 0.0537, PSNR (dB) = 12.70\n",
      "MSE = 0.0557, PSNR (dB) = 12.54\n",
      "MSE = 0.0180, PSNR (dB) = 17.45\n",
      "MSE = 0.0571, PSNR (dB) = 12.44\n",
      "MSE = 0.0615, PSNR (dB) = 12.11\n",
      "MSE = 0.0780, PSNR (dB) = 11.08\n",
      "MSE = 0.0411, PSNR (dB) = 13.86\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": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3857\n",
      "0.2428\n",
      "0.3230\n",
      "0.4340\n",
      "0.1804\n",
      "0.4266\n",
      "0.1498\n",
      "0.0997\n",
      "0.1599\n",
      "0.5370\n",
      "Wall time: 213 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/imamun93/Image-Similarities-using-SSIM/blob/master/blog5_ssim.ipynb\n",
    "import numpy as np\n",
    "from tensorflow.image import ssim\n",
    "from tensorflow import convert_to_tensor\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = np.asarray(ssim(convert_to_tensor(original), convert_to_tensor(contrast), max_val=1)) # <class 'tensorflow.python.framework.ops.EagerTensor'>\n",
    "    print(f'{a[0]:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6480\n",
      "0.4463\n",
      "0.8745\n",
      "0.6208\n",
      "0.1592\n",
      "0.9007\n",
      "0.4673\n",
      "0.2496\n",
      "0.1482\n",
      "0.8138\n",
      "Wall time: 165 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# https://github.com/tensorflow/models/blob/master/research/compression/image_encoder/msssim.py\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from scipy.ndimage.filters import convolve\n",
    "\n",
    "def _FSpecialGauss(size, sigma):\n",
    "  \"\"\"Function to mimic the 'fspecial' gaussian MATLAB function.\"\"\"\n",
    "  radius = size // 2\n",
    "  offset = 0.0\n",
    "  start, stop = -radius, radius + 1\n",
    "  if size % 2 == 0:\n",
    "    offset = 0.5\n",
    "    stop -= 1\n",
    "  x, y = np.mgrid[offset + start:stop, offset + start:stop]\n",
    "  assert len(x) == size\n",
    "  g = np.exp(-((x**2 + y**2)/(2.0 * sigma**2)))\n",
    "  return g / g.sum()\n",
    "\n",
    "\n",
    "def _SSIMForMultiScale(img1, img2, max_val=255, filter_size=11,\n",
    "                       filter_sigma=1.5, k1=0.01, k2=0.03):\n",
    "  \"\"\"Return the Structural Similarity Map between `img1` and `img2`.\n",
    "  This function attempts to match the functionality of ssim_index_new.m by\n",
    "  Zhou Wang: http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "  Returns:\n",
    "    Pair containing the mean SSIM and contrast sensitivity between `img1` and\n",
    "    `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  img1 = img1.astype(np.float64)\n",
    "  img2 = img2.astype(np.float64)\n",
    "  _, height, width, _ = img1.shape\n",
    "\n",
    "  # Filter size can't be larger than height or width of images.\n",
    "  size = min(filter_size, height, width)\n",
    "\n",
    "  # Scale down sigma if a smaller filter size is used.\n",
    "  sigma = size * filter_sigma / filter_size if filter_size else 0\n",
    "\n",
    "  if filter_size:\n",
    "    window = np.reshape(_FSpecialGauss(size, sigma), (1, size, size, 1))\n",
    "    mu1 = signal.fftconvolve(img1, window, mode='valid')\n",
    "    mu2 = signal.fftconvolve(img2, window, mode='valid')\n",
    "    sigma11 = signal.fftconvolve(img1 * img1, window, mode='valid')\n",
    "    sigma22 = signal.fftconvolve(img2 * img2, window, mode='valid')\n",
    "    sigma12 = signal.fftconvolve(img1 * img2, window, mode='valid')\n",
    "  else:\n",
    "    # Empty blur kernel so no need to convolve.\n",
    "    mu1, mu2 = img1, img2\n",
    "    sigma11 = img1 * img1\n",
    "    sigma22 = img2 * img2\n",
    "    sigma12 = img1 * img2\n",
    "\n",
    "  mu11 = mu1 * mu1\n",
    "  mu22 = mu2 * mu2\n",
    "  mu12 = mu1 * mu2\n",
    "  sigma11 -= mu11\n",
    "  sigma22 -= mu22\n",
    "  sigma12 -= mu12\n",
    "\n",
    "  # Calculate intermediate values used by both ssim and cs_map.\n",
    "  c1 = (k1 * max_val) ** 2\n",
    "  c2 = (k2 * max_val) ** 2\n",
    "  v1 = 2.0 * sigma12 + c2\n",
    "  v2 = sigma11 + sigma22 + c2\n",
    "  ssim = np.mean((((2.0 * mu12 + c1) * v1) / ((mu11 + mu22 + c1) * v2)))\n",
    "  cs = np.mean(v1 / v2)\n",
    "  return ssim, cs\n",
    "\n",
    "\n",
    "def MultiScaleSSIM(img1, img2, max_val=255, filter_size=11, filter_sigma=1.5,\n",
    "                   k1=0.01, k2=0.03, weights=None):\n",
    "  \"\"\"Return the MS-SSIM score between `img1` and `img2`.\n",
    "  This function implements Multi-Scale Structural Similarity (MS-SSIM) Image\n",
    "  Quality Assessment according to Zhou Wang's paper, \"Multi-scale structural\n",
    "  similarity for image quality assessment\" (2003).\n",
    "  Link: https://ece.uwaterloo.ca/~z70wang/publications/msssim.pdf\n",
    "  Author's MATLAB implementation:\n",
    "  http://www.cns.nyu.edu/~lcv/ssim/msssim.zip\n",
    "  Arguments:\n",
    "    img1: Numpy array holding the first RGB image batch.\n",
    "    img2: Numpy array holding the second RGB image batch.\n",
    "    max_val: the dynamic range of the images (i.e., the difference between the\n",
    "      maximum the and minimum allowed values).\n",
    "    filter_size: Size of blur kernel to use (will be reduced for small images).\n",
    "    filter_sigma: Standard deviation for Gaussian blur kernel (will be reduced\n",
    "      for small images).\n",
    "    k1: Constant used to maintain stability in the SSIM calculation (0.01 in\n",
    "      the original paper).\n",
    "    k2: Constant used to maintain stability in the SSIM calculation (0.03 in\n",
    "      the original paper).\n",
    "    weights: List of weights for each level; if none, use five levels and the\n",
    "      weights from the original paper.\n",
    "  Returns:\n",
    "    MS-SSIM score between `img1` and `img2`.\n",
    "  Raises:\n",
    "    RuntimeError: If input images don't have the same shape or don't have four\n",
    "      dimensions: [batch_size, height, width, depth].\n",
    "  \"\"\"\n",
    "  if img1.shape != img2.shape:\n",
    "    raise RuntimeError('Input images must have the same shape (%s vs. %s).',\n",
    "                       img1.shape, img2.shape)\n",
    "  if img1.ndim != 4:\n",
    "    raise RuntimeError('Input images must have four dimensions, not %d',\n",
    "                       img1.ndim)\n",
    "\n",
    "  # Note: default weights don't sum to 1.0 but do match the paper / matlab code.\n",
    "  weights = np.array(weights if weights else\n",
    "                     [0.0448, 0.2856, 0.3001, 0.2363, 0.1333])\n",
    "  levels = weights.size\n",
    "  downsample_filter = np.ones((1, 2, 2, 1)) / 4.0\n",
    "  im1, im2 = [x.astype(np.float64) for x in [img1, img2]]\n",
    "  mssim = np.array([])\n",
    "  mcs = np.array([])\n",
    "  for _ in range(levels):\n",
    "    ssim, cs = _SSIMForMultiScale(\n",
    "        im1, im2, max_val=max_val, filter_size=filter_size,\n",
    "        filter_sigma=filter_sigma, k1=k1, k2=k2)\n",
    "    mssim = np.append(mssim, ssim)\n",
    "    mcs = np.append(mcs, cs)\n",
    "    filtered = [convolve(im, downsample_filter, mode='reflect')\n",
    "                for im in [im1, im2]]\n",
    "    im1, im2 = [x[:, ::2, ::2, :] for x in filtered]\n",
    "  return (np.prod(mcs[0:levels-1] ** weights[0:levels-1]) *\n",
    "          (mssim[levels-1] ** weights[levels-1]))\n",
    "\n",
    "n = 10\n",
    "for i in range(n):\n",
    "    original = x_test[i].reshape(1, 28, 28, 1)\n",
    "    contrast = decoded_imgs[i].reshape(1, 28, 28, 1)\n",
    "    a = MultiScaleSSIM(original, contrast, max_val=1.0)\n",
    "    print(f'{a:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Simple Autoencoders in Keras - 784-2-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
}