Deep Autoencoders in Keras - 784-256-128-64-32-16-32-64-128-256-784-tf1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import os\n",
    "#os.environ['KERAS_BACKEND'] = 'plaidml.keras.backend' # https://github.com/plaidml/plaidml/blob/master/docs/index.md\n",
    "\n",
    "from tensorflow.keras.models import Model\n",
    "from tensorflow.keras.layers import Input, Dense\n",
    "from tensorflow.keras.datasets import mnist\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# keras version still end up import tf.keras version which is incompatible, import directly instead\n",
    "from tensorflow.keras.callbacks import TensorBoard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "encoding_dim = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_img = Input(shape=(784,))\n",
    "encoded = Dense(256, activation='relu')(input_img) # apply to previous layer = input_img\n",
    "encoded = Dense(128, activation='relu')(encoded)\n",
    "encoded = Dense(64, activation='relu')(encoded)\n",
    "encoded = Dense(32, activation='relu')(encoded)\n",
    "encoded = Dense(encoding_dim, activation='relu')(encoded) # the bottleneck = encoding_dim\n",
    "\n",
    "decoded = Dense(32, activation='relu')(encoded)\n",
    "decoded = Dense(64, activation='relu')(decoded)\n",
    "decoded = Dense(128, activation='relu')(decoded)\n",
    "decoded = Dense(256, activation='relu')(decoded)\n",
    "decoded = Dense(784, activation='sigmoid')(decoded)\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[-5](encoded_input)\n",
    "decoder_layer = autoencoder.layers[-4](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-3](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-2](decoder_layer)\n",
    "decoder_layer = autoencoder.layers[-1](decoder_layer) # extract last layer\n",
    "\n",
    "decoder = Model(encoded_input, decoder_layer)"
   ]
  },
  {
   "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": [
      "1595257932\n",
      "Epoch 1/50\n",
      "  2/235 [..............................] - ETA: 3:50 - loss: 0.6930WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.997496). Check your callbacks.\n",
      "  3/235 [..............................] - ETA: 2:38 - loss: 0.6925WARNING:tensorflow:Method (on_train_batch_end) is slow compared to the batch update (0.109996). Check your callbacks.\n",
      "235/235 [==============================] - 66s 281ms/step - loss: 0.2632 - val_loss: 0.1898\n",
      "Epoch 2/50\n",
      "235/235 [==============================] - 66s 281ms/step - loss: 0.1725 - val_loss: 0.1601\n",
      "Epoch 3/50\n",
      "235/235 [==============================] - 66s 282ms/step - loss: 0.1531 - val_loss: 0.1460\n",
      "Epoch 4/50\n",
      "235/235 [==============================] - 71s 301ms/step - loss: 0.1430 - val_loss: 0.1389\n",
      "Epoch 5/50\n",
      "235/235 [==============================] - 74s 314ms/step - loss: 0.1373 - val_loss: 0.1337\n",
      "Epoch 6/50\n",
      "235/235 [==============================] - 75s 319ms/step - loss: 0.1321 - val_loss: 0.1293\n",
      "Epoch 7/50\n",
      "235/235 [==============================] - 80s 342ms/step - loss: 0.1287 - val_loss: 0.1272\n",
      "Epoch 8/50\n",
      "235/235 [==============================] - 85s 364ms/step - loss: 0.1252 - val_loss: 0.1227\n",
      "Epoch 9/50\n",
      "235/235 [==============================] - 92s 390ms/step - loss: 0.1223 - val_loss: 0.1206\n",
      "Epoch 10/50\n",
      "235/235 [==============================] - 90s 384ms/step - loss: 0.1204 - val_loss: 0.1187\n",
      "Epoch 11/50\n",
      "235/235 [==============================] - 90s 382ms/step - loss: 0.1187 - val_loss: 0.1177\n",
      "Epoch 12/50\n",
      "235/235 [==============================] - 88s 376ms/step - loss: 0.1172 - val_loss: 0.1158\n",
      "Epoch 13/50\n",
      "235/235 [==============================] - 90s 385ms/step - loss: 0.1156 - val_loss: 0.1144\n",
      "Epoch 14/50\n",
      "235/235 [==============================] - 90s 381ms/step - loss: 0.1142 - val_loss: 0.1137\n",
      "Epoch 15/50\n",
      "235/235 [==============================] - 80s 341ms/step - loss: 0.1130 - val_loss: 0.1124\n",
      "Epoch 16/50\n",
      "235/235 [==============================] - 73s 311ms/step - loss: 0.1121 - val_loss: 0.1113\n",
      "Epoch 17/50\n",
      "235/235 [==============================] - 69s 293ms/step - loss: 0.1111 - val_loss: 0.1107\n",
      "Epoch 18/50\n",
      "235/235 [==============================] - 65s 277ms/step - loss: 0.1103 - val_loss: 0.1098\n",
      "Epoch 19/50\n",
      "235/235 [==============================] - 66s 282ms/step - loss: 0.1093 - val_loss: 0.1098\n",
      "Epoch 20/50\n",
      "235/235 [==============================] - 60s 255ms/step - loss: 0.1082 - val_loss: 0.1089\n",
      "Epoch 21/50\n",
      "235/235 [==============================] - 51s 218ms/step - loss: 0.1073 - val_loss: 0.1069\n",
      "Epoch 22/50\n",
      "235/235 [==============================] - 64s 271ms/step - loss: 0.1048 - val_loss: 0.1052\n",
      "Epoch 26/50\n",
      "235/235 [==============================] - 70s 296ms/step - loss: 0.1044 - val_loss: 0.1047\n",
      "Epoch 27/50\n",
      "235/235 [==============================] - 80s 339ms/step - loss: 0.1038 - val_loss: 0.1049\n",
      "Epoch 28/50\n",
      "235/235 [==============================] - 79s 338ms/step - loss: 0.1035 - val_loss: 0.1048\n",
      "Epoch 29/50\n",
      "235/235 [==============================] - 49s 211ms/step - loss: 0.1032 - val_loss: 0.1035\n",
      "Epoch 30/50\n",
      "235/235 [==============================] - 64s 273ms/step - loss: 0.1028 - val_loss: 0.1033\n",
      "Epoch 31/50\n",
      "235/235 [==============================] - 62s 265ms/step - loss: 0.1025 - val_loss: 0.1031\n",
      "Epoch 32/50\n",
      "235/235 [==============================] - 57s 242ms/step - loss: 0.1022 - val_loss: 0.1029\n",
      "Epoch 33/50\n",
      "235/235 [==============================] - 51s 215ms/step - loss: 0.1020 - val_loss: 0.1026\n",
      "Epoch 34/50\n",
      "235/235 [==============================] - 59s 252ms/step - loss: 0.1017 - val_loss: 0.1029\n",
      "Epoch 35/50\n",
      "235/235 [==============================] - 57s 242ms/step - loss: 0.1015 - val_loss: 0.1026\n",
      "Epoch 36/50\n",
      "235/235 [==============================] - 53s 226ms/step - loss: 0.1011 - val_loss: 0.1022\n",
      "Epoch 37/50\n",
      "235/235 [==============================] - 61s 262ms/step - loss: 0.1010 - val_loss: 0.1017\n",
      "Epoch 38/50\n",
      "235/235 [==============================] - 57s 241ms/step - loss: 0.1007 - val_loss: 0.1025\n",
      "Epoch 39/50\n",
      "235/235 [==============================] - 58s 248ms/step - loss: 0.1006 - val_loss: 0.1018\n",
      "Epoch 40/50\n",
      "235/235 [==============================] - 58s 249ms/step - loss: 0.1003 - val_loss: 0.1017\n",
      "Epoch 41/50\n",
      "235/235 [==============================] - 47s 199ms/step - loss: 0.1001 - val_loss: 0.1014\n",
      "Epoch 42/50\n",
      "235/235 [==============================] - 57s 241ms/step - loss: 0.0999 - val_loss: 0.1012\n",
      "Epoch 43/50\n",
      "235/235 [==============================] - 59s 251ms/step - loss: 0.0998 - val_loss: 0.1012\n",
      "Epoch 44/50\n",
      "235/235 [==============================] - 43s 184ms/step - loss: 0.0996 - val_loss: 0.1009\n",
      "Epoch 45/50\n",
      "235/235 [==============================] - 39s 166ms/step - loss: 0.0995 - val_loss: 0.1007\n",
      "Epoch 46/50\n",
      "235/235 [==============================] - 41s 173ms/step - loss: 0.0993 - val_loss: 0.1014\n",
      "Epoch 47/50\n",
      "235/235 [==============================] - 39s 166ms/step - loss: 0.0993 - val_loss: 0.1007\n",
      "Epoch 48/50\n",
      "235/235 [==============================] - 44s 187ms/step - loss: 0.0990 - val_loss: 0.1009\n",
      "Epoch 49/50\n",
      "235/235 [==============================] - 35s 148ms/step - loss: 0.0989 - val_loss: 0.1007\n",
      "Epoch 50/50\n",
      "235/235 [==============================] - 40s 171ms/step - loss: 0.0988 - val_loss: 0.1004\n",
      "1595257932\n",
      "Wall time: 53min 24s\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": 7,
   "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+WH4yJAAAIvUlEQVR4nO3bS4iWdRvH8WtERuZQjg6TM20saZkgoiBUoh3AQSOZkA5UgogEIrYIElxpEtjCkBKJRB2sRW2jbUsx9V1oKJopDHigHF0EqQvjbvuiIzQv95XX4/v5LOXhy+38n/8cfsx0NU0TAAAAANQy42E/AAAAAAD3M9oAAAAAFGS0AQAAACjIaAMAAABQkNEGAAAAoCCjDQAAAEBBM6fz4oGBgWZkZCTrWWLmzGk9zrRcunQprR0RcevWrdR+0zRdbXS6u7ub3t7eNlJTeuKJJ9LaM2bkboznz59P7UfEZNM0Q22Euru7m56enjZSU5o3b15a+48//khrR0T89ttvqf227mJPT08ze/bsNlJTGhgYSGvfvHkzrR0Rcf369dR+tHgXH3/88SbzvmS+R86dO5fWjoj4888/U/tt3cVZs2alfl3MPMPbt2+ntSMifv/999R+tHgXH3vssWZoqJXUlDLP8cKFC2ntiM65iwMDA83w8HAbqSn19/entScnJ9PaERETExOp/WjxLvb19TVz5sxpIzWlwcHBtPa1a9fS2hH539+0dRf7+/ubzI9zphs3bqT2sz+fxgPu4rRWkpGRkTh8+HBrT3SvzDfHW2+9ldaOiDh58mRqvy29vb2xfPnytP7mzZvT2n19fWntiIgXXnghtR8RrX3F7enpieeee66t3H22bt2a1v7xxx/T2hERn376aWq/LbNnz4533nknrT82NpbW/vrrr9PaERH79+9P7UeLd3HevHmxZ8+etnL3efXVV9PamZ9DIiKOHj2a2m9Lb29vvPTSS2n90dHRtPaZM2fS2hERn332WWo/WryLQ0ND8cknn7SVu8+qVavS2pnvkYiIY8eOpfbbMjw8HAcOHEjrP//882ntr776Kq0dEbFp06bUfrR4F+fMmZP688CGDRvS2rt27UprR0R88cUXqf22DA4Oxvbt29P6f/31V1o7+3vUf+F7mynvoj+PAgAAACjIaAMAAABQkNEGAAAAoCCjDQAAAEBBRhsAAACAgow2AAAAAAUZbQAAAAAKMtoAAAAAFGS0AQAAACjIaAMAAABQkNEGAAAAoCCjDQAAAEBBRhsAAACAgow2AAAAAAUZbQAAAAAKMtoAAAAAFGS0AQAAACjIaAMAAABQkNEGAAAAoCCjDQAAAEBBRhsAAACAgrqapvnnL+7q+ucv/h9M51mmq6urK60dEXHlypW09ujoaJw6daqV/4AzfLA33ngjtf/tt9/+p2maJW20nOOD/fLLL2ntsbGx+Pnnn93F5DNcv359an98fLxj7uLVq1fT2k8++WRaOyLiu+++S2tv27YtLl686C4m38VLly6l9hcsWNAxd7GTz/GHH35Ia2/dujUuXLjQEXfx2rVrae2RkZG0dkTEvn37UvubN292FyP/Lp4+fTqt/eabb8aZM2c64i528hkeO3Ystb9s2bIp76LftAEAAAAoyGgDAAAAUJDRBgAAAKAgow0AAABAQUYbAAAAgIKMNgAAAAAFGW0AAAAACjLaAAAAABRktAEAAAAoyGgDAAAAUJDRBgAAAKAgow0AAABAQUYbAAAAgIKMNgAAAAAFGW0AAAAACjLaAAAAABRktAEAAAAoyGgDAAAAUJDRBgAAAKAgow0AAABAQUYbAAAAgIKMNgAAAAAFzXzYD/Df9u7dm9Zet25dWjsiYtu2bWnty5cvp7XbtnDhwrT2hx9+mNaOiDhx4kRqv5O8+OKLae2DBw+mtSMixsbG0tq//vprWrttX375ZVr75ZdfTmtHRKxcuTK1Pz4+ntpv00cffZTWfv3119PaERE3btxIa9+9ezet3bbR0dG09pEjR9LaERFr1qxJ7XeSFStWpLX7+/vT2hERBw4cSGtPTk6mtdu2Y8eOtPahQ4fS2hERExMTqf1Okvmzxtq1a9PaERHnz59Pa9+5cyet3bbMr4uvvPJKWjsi4ptvvkntP4jftAEAAAAoyGgDAAAAUJDRBgAAAKAgow0AAABAQUYbAAAAgIKMNgAAAAAFGW0AAAAACjLaAAAAABRktAEAAAAoyGgDAAAAUJDRBgAAAKAgow0AAABAQUYbAAAAgIKMNgAAAAAFGW0AAAAACjLaAAAAABRktAEAAAAoyGgDAAAAUJDRBgAAAKAgow0AAABAQUYbAAAAgIKMNgAAAAAFGW0AAAAACpo5nRf39fXFokWLsp4lPvjgg7T24sWL09oREXv27Elr//TTT6215s+fH9u3b2+td69NmzaltZctW5bWjoiYMaNzNszh4eFYv359Wn/37t1p7du3b6e1IyIWLlyY1r569WprrcHBwVizZk1rvXu9//77ae133303rR0R8eyzz6b22/T000/Hzp070/qZH+stW7aktSMiJicn09p3795trTUyMhIbN25srXevjz/+OK19+fLltHZExOjoaGr/7NmzrbWeeuqp1Lv43nvvpbWXLFmS1o6IWLp0aVr7+PHjrbU6+XvUO3fupLUjIt5+++3UfpueeeaZ2Lt3b1p/9erVae02389T+f7771P7bZk7d27qx/nIkSNp7ddeey2tHRGxcuXK1P7nn38+5b93zk+pAAAAAP9HjDYAAAAABRltAAAAAAoy2gAAAAAUZLQBAAAAKMhoAwAAAFCQ0QYAAACgIKMNAAAAQEFGGwAAAICCjDYAAAAABRltAAAAAAoy2gAAAAAUZLQBAAAAKMhoAwAAAFCQ0QYAAACgIKMNAAAAQEFGGwAAAICCjDYAAAAABRltAAAAAAoy2gAAAAAUZLQBAAAAKMhoAwAAAFCQ0QYAAACgoK6maf75i7u6rkfERN7j8ADzm6YZaiPkDB8q59j5nOGjwTl2Pmf4aHCOnc8ZPhqcY+dzho+GKc9xWqMNAAAAAP8Ofx4FAAAAUJDRBgAAAKAgow0AAABAQUYbAAAAgIKMNgAAAAAFGW0AAAAACjLaAAAAABRktAEAAAAoyGgDAAAAUNDfb4bLJEmUYWkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x576 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAABwCAYAAACkaY2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debTVZfXH8U3lUKgMggiIgqCIIBAqFsrgyjTFHFIrzXJCl+UKnMhiSblyqWnmsFYtWi6zhFKXA1oOGIo5kIiCoOCEyCwyXAecUBv4/eFq/z57c+4RLufe+73nvF9/7ePzcO7X8z3P9/uc73r2s1tt2LDBAAAAAAAAUCyfa+4DAAAAAAAAwMZ4aAMAAAAAAFBAPLQBAAAAAAAoIB7aAAAAAAAAFBAPbQAAAAAAAAqIhzYAAAAAAAAF9IXN6dyqVSvqgzeTDRs2tKrE+3AOm1Xdhg0bOlbijTiPzYexWBUYi1WAsVgVGItVgLFYFRiLVYCxWBVKjkVW2gBNZ2lzHwAAM2MsAkXBWASKgbEIFEPJschDGwAAAAAAgALioQ0AAAAAAEAB8dAGAAAAAACggHhoAwAAAAAAUECbVT0KAFA7WrX6/yIEGzZQSAAAAABoaqy0AQAAAAAAKCAe2gAAAAAAABQQ6VGoKE2n2HbbbUNbhw4dPO7Tp4/He+yxR+jXt29fjzt37uzx5z//+dDv/vvv9/iee+4JbXV1dR7/5z//2aRjR/GQntP49DNu165daBswYIDHXbt29fjll18O/RYsWODxBx98ENoYfwAAAGbbb799eK1z23//+98ef/TRR012TGgZWGkDAAAAAABQQDy0AQAAAAAAKCAe2gAAAAAAABQQe9qgolq3bu3xiSeeGNpGjx7tcfv27T3eZpttQr/Pfe5zJdu+8IX4dR06dKjHZ511Vmg74YQTPF68eLHH7ItSPHq+L7zwwtB27LHHenzjjTd6fNNNN4V+nNeG072idtppp9A2atQoj/fbbz+P81jUPW4mTJgQ2h588EGPNV8bjUfHlFnct+i///2vxw0dN/p+ea8x/dv5fOvfRuPL+8p17NjR49WrV3v8ySefNNkxoenoOM24ZzZcuc+1vn5c+6qTnmONv/jFL4Z+3bt39zjPc3X+dMcdd3i8bt260E/3C2S/m9rEShsAAAAAAIAC4qENAAAAAABAAZEehS2S0yROPvlkj3/605+GNi0nrP8uL/P717/+5bEuB9xhhx1CPy2bt/vuu4e2wYMHe7xkyRKPWRJcPJpeoSWmzWKZ+B133LHJjqmWaApi/vx1Sa/KZbx1bA4aNCi0Pfnkkx6/9dZbDT1MlKBjR8dK586dQ7+6ujqPV61a5fHmpKvp0u8vfelL9f4tvbZraqqZ2ccff7zJfw8No+lpY8eODW2nn366x/fdd5/Hmrpsxn2yPjoGtt5669Cm9ydte/vtt0M/ndOUG3/lUnC22morj7t06eLxl7/85dBPr8uPPfZYaNP0OJ2DVfu512tm/n+tL+0zn2u9xul8Nc+H9XpHOkux6bnX8WUWf1/069cvtA0fPtxjvS/m70yvXr081t9CZmY9evTw+M033/R4+vTpod8777xT//8AagIrbQAAAAAAAAqIhzYAAAAAAAAF1CTpUbrsTJcj5uWf9fXLqn35ZkuSK4csWrTI45deeim07bzzzh7r0tx58+aFfppOseeee3qcK0Tp7uy5SsZhhx3m8Z133ukxO/gXj1Y0yek4uqT4gQce8JhrwJbRa61WctOqa2ZmPXv2LPlv8rjXZfht27YNbR9++KHHWllK0wTQMG3atPH4pJNO8jinqGlFiilTpjTob+mY01iXfZuZvf/++x4vXbq0QX8LDadL9E855ZTQ1rVrV481hThXG8vpj/iUppL2798/tPXt29fj7bbbzuMZM2aEfq+++qrHOlbKzU3y+dHUiz322MPjI488st73yKmKWplG03iq/d6qn3P+DaKfq1ZSzGlnel/Uf5MroT7yyCMe51QX5qJNL59vPV+63YJWyjQzO/vssz3WVCaz+NtDU+Xyd0GvCTonMovfoW7dunm8du3a0E+vA1yjaxMrbQAAAAAAAAqIhzYAAAAAAAAFxEMbAAAAAACAAmryPW00T3To0KGhn+aN6v4nZjG3b86cOR7PnTs39Fu5cqXHmqeb80c1bzfnC2uOou73kHMU9Zg0P7jU36tWuWTl7NmzPR4/fnxo089V86vffffd0O+TTz7x+PHHH/dY928wi3vc5HOopYXLlc5E89P9UPK5mj9/vscLFixosmOqdjpehg0b5vHAgQNDP903Sq+ZeU+b1q1be5zH6UUXXeSxXuN//OMfh36UA/9s+Tp3wAEHeDxq1CiPc868jh3Nu28oLYm6//77h7Zly5Z5/Pzzz4e29evXb/HfRnlHH320x7vuumto03Grpd9rZb6ypXRfihEjRoQ23X9P71u5TK+OzU3dlyL30zGs1+g8b9Y9BnMZ43J7u9SKfB/T3ydHHHGEx6eddlrop5+z7oWSS37r3or6fmYbl4JH49DfHXmvGr1W6t5QumeRmVmHDh08XrJkSWhbsWKFx/p90r3FzOLvnLxXjf5ufeGFFzzOZeJrdR8bHVedOnUKbbqXmM6H9t1339BPz02+Jq9Zs8ZjPZ/5PGk59jy30b56bc3nbEv3DGOlDQAAAAAAQAHx0AYAAAAAAKCAmiQ9SpcDaTpN586dQz8tN6tLDs3iElBdGlpXVxf66fvrEqWcxqOpU3nJuf4tXeKWS9nq0lMtt2pm9tprr1ktyEu/NE1M09jM6i/pnpdmaz9dxqalMs3icjctnWlmNnXqVI/zuUex6HUgj7HrrrvOY02bw5Zp166dx5pWo8uAzWIpyvpKPpvFa2hecq5/69hjjy35b8xieWLOdWk5xUFLkeqy4d/+9reh38KFCyt6HJoCl9Oj9DvDtbfx5bTtiy++2OP8fdG505VXXulxtZd5bqh8jdL00fy919Sp++67z+PXX3899KtEeqIel6Z15NQBTa/IKfyaqlir5z/fq3r37u2xpjPpuTWL1z9NDc7vpykaBx54YGi7//77Pa7Vz79ScnqfpirqdU7TZ8zi9VHT1TR11Cyeq1mzZoW2N954o+Rx5DmMpn9rmo1ZvCboPVN/p5rV1vdEx9jll1/u8THHHBP66fwy3++U/s784IMPQpteG9977z2PdeuG/P45TW7ixIkeT5s2zeN8/ddrckPOJyttAAAAAAAACoiHNgAAAAAAAAXEQxsAAAAAAIACavI9bTSv76677gr9NDf7kEMOCW1aKlZzBXOuqe6Fo7mBObdeyy6uXr06tGn5r+7du3usuXNmZv369fM472lz2WWXeVxL5TQ3dQ8Dzf3M+aj6umPHjh6PHTs29NNzr2VmzWLp8VrKA20J8j4Bus9JLpM4d+7cJjmmapdLker1atCgQR7rfl5mcc8qve4uX7489NNc3/bt24c2zU3W/U6OPPLI0O/cc8/1+Oqrrw5ttXQNLecrX/lKeK33Sc3FnjRpUuhX6b1ltJR0//79Q5veP2u1lHBT0jmKmVmPHj3q7au5+7lkKTaW90g49NBDPdb5n1ncE0P3O8j7UjREvmfutttuHh9++OEea8lqM7NXXnml5PGZxWt7rc6Rcol0vQfts88+HufPZ+nSpSXb8ljU/W6GDh0a2v7+9797XIl9jmqNjonjjjsutF111VUe5336lJZ3/vOf/+xxvn/q3pp5rxrdx0jvd/meW2686fym3N6B1Ux/65mZ3XvvvR7rHDXvG6WfkZ6bfJ50Lxk97/m1jtm8R5g+X8j73Rx11FEe636ueU+bLcVKGwAAAAAAgALioQ0AAAAAAEABNUl6lNKlTCtXrgxtWt73lltuCW26HFuX8Oelb7oESpdR6dJxs7hMWJe+ZT/72c9KHoNZXNJYrtQYPqXLGXUZYV7upukUY8aM8bhbt271vvfkyZPD61xSD8WRx8rRRx/tcV5SquUU0XBaAtPM7KKLLvJY0wzz569pUA8//LDHd999d+inY3jYsGGh7bvf/a7HXbp08TinYv3kJz/x+IEHHght8+fPt1ql4+W8884LbXpNffDBBz3O5SgrQf+WpjTmJcSaRsyy/8ah988hQ4aENk1BzMvrddm2poijtLZt24bXmn6UU/P1s9WSwQ1NcdBzvO2224a20047zeM+ffp4nFOxHn30UY/r6upCW30pGdVOP9eePXuGNp1j6rXriSeeCP2mTp3qsV4XTz311NBv8ODBHuf7oqZXMF/9bDlFUMfApZdeGto0xWXt2rUe629MM7Obb77ZY/2NmMdDuTRfvT/rmCqXklxL460cPU9/+tOfQtvAgQM91vT+XK5br3H63EDTtM3ivES3YjGLzxFGjBjhcb7u6vcgbzmg/y/lvgdbeu5ZaQMAAAAAAFBAPLQBAAAAAAAooCZPjypHlxHl1ClNk8jL5FRDlnzmpW+6dEqXGufdqPWYdOfxfBzYWLmlv8cff7zHJ5xwgsf5M33yySc9vvLKK0Mbyw+Lq1evXuG17hp/6623hjaW8DecVuL65S9/GdryTv3/k9NIf/Ob33h82223eZzPi45nHZdmMRX1ggsu8DinHmhapFbfMzM75phjPK61sb333nt7nCs16Xm49tprPdZKFZWiy8C//e1ve5yv35qa1RjHgThHOf3000ObjkVNFzcz++Mf/+gxc5TS9PPLKdlaHShfh7RKiH62eX5Z3/Ur99Pl95pWamb2ta99zWP9LsyaNSv00wpF77///iYdR7UrV7lUU8gWLVrk8cSJE0O/l156yWM9T7mimKauafqEWawAplV1a/W8/I+eE61EedZZZ4V+OpfIn9nf/vY3jy+++GKPc5XZhvxGzGkx+htR78e1fh5LyVth6DziwAMPDG0639CUqNGjR4d++ptBUxpz+r3+7shV9jQVS6vG5XOtzx7yHFgr3ZYbz/pdash3hJU2AAAAAAAABcRDGwAAAAAAgALioQ0AAAAAAEABFWpPm3I096vSefI5r0z3gtA9BPKeNnfeeafHixcvrugxVSP9nDW3UferMDO78MILPdZSbLk0+7hx4zzOufsoFs0FHTt2bGjTHE8dU2bsibElhg8f7vHBBx8c2jRXV8eOljI1M5s0aZLHm7q/UL5OTpgwwWPN4z/llFNCP81B3n///UOb7ttQ7WM951GPGjXKY70emsUymY1dFl1LZmq52vy9eOyxxzxm35TGoTn5mo+fvfvuu+G17vWA0vR+tMsuu4Q23Wch35t0TtOjRw+Pc2nw+q5feQ8GfZ3nSHocWi76rrvuCv1030XupZ/SeWj+zPWz1OvpihUrQj+9x+m1MF+fy+1xot+tefPmeaz7ctSirl27enzmmWd6PGbMmNBP5wTTp08Pbeeff77HWvK7oXT+mvdU1b1Yud+Vl+c2I0eO9DiPRR0Hum9NvsZ9/PHHJd8/7y+lpbwPOOCA0Lbvvvt6rGM473ml4/75558PbVOmTPG4vv3NzCj5DQAAAAAAUJV4aAMAAAAAAFBALSY9qjHlMmTf+c53PNYSj3mJ5A033OBxrS9p3BS6dG2vvfby+KSTTgr9dNmoLunN5SyfffbZSh8iGokuEde0HTOzVatWeZzLRWPT6XJhM7MzzjjD4+233z606RLNOXPmeHzOOeeEfpUoua7lGjX9Sss9msXS0XmprJZLrfb0qB133DG81mW9OcVB06Mqnf6QlwbvvPPOHuuS31xGdfbs2R5T9rRxaPpNXnKu50ZLPptVZjxXO/3e5znfggULPM5pEnvvvbfH119/vcedOnUK/fTf6djJKfaaqjN48OB6j1GP6aGHHgr9NHUDnyqX/qbnStOj8v1TU0U11lQNs3gO89gbMmSIxzNnzvQ4p/NUe8pNHkdXXHGFx8cdd5zHeX6jJewff/zx0LZu3bpKHmI4Bzn9W+9x3O/Ky3OK9957z+P82elnvuuuu3qcfz/o7/fvfe97HmvKk1nc9kRTGs1i2XY9xjye9X6gKVtmcd6jc7FKj19W2gAAAAAAABQQD20AAAAAAAAKiPQoM+vTp094fcEFF3isS/Gvvfba0G/lypWNe2BVRpc37rPPPh7rUm+zuExOq1/8/ve/D/3yMkUUl+7WrhWEzMymTZvmsS55xebJ42jQoEEe5yXIulT70ksv9ThXaKsEXb6qS0/zMam8LDUvq61m+Tx269bN43zN089F40os085LiDWNdenSpR4/8cQToR9juHHoeDnooIM81qXdZjFV+5prrgltLN//bLqc/YUXXghtEydO9Hi//fYLbb169fJYK+C0bds29NN0Nr0OL1y4MPTbc889PdZ0crOYcqrX75xaw/nemN6PunTpEtr0Wjts2DCPtVqbWfyc9XdAud8E+XeGpm8cccQRHucKb2+//bbH1Zgqle/tOufXsZLnC5qCotXUzGIFoLfeeqvkvyn1+n/yNVWPI6cc6j2Z8VZe/v5qaqemSpnFlPgDDzzQY00rNIup9PpdyudJ74v5/Oo1QY9R0xvNzC6//HKPJ0+eHNrKVYqrJFbaAAAAAAAAFBAPbQAAAAAAAAqIhzYAAAAAAAAFVLN72mh52V//+tehTcsTP/bYYx7fcsstoV+lS6xWm5yr2qZNG48HDhzocc5H1bxgLZ359NNPh37N9fnn/y/yWEvT8zp+/HiP8+d38803e8yYarjOnTuH15qHnT9XHWNPPfWUx43xXda9rHRvI70Gm8Vc4jVr1oS2Whpj+f9Vc6W1bKWZ2amnnuqx7iWzfPny0E/3CNLxl0ura2nNs846K7R961vfKvl+zzzzTOjHXmONQ/dGOfnkkz3Oez3ofhtLlixp9OOqNjr+1q9fH9p0/6YZM2aENr3f7bjjjh7nPar0PTVu37596Pfzn//cY52TmpmtWrXK4/pKzaI0HS/du3cPbVraW/eE071uzMyeffbZku+Xx6Le43T+axa/LyNGjPA4z3P1up73eqtGuveIlu7WPU7M4v57eX+pSy65pOS/y/cm3bNE99LRPU7M4txkzpw5oe3RRx/1WO/VeU8VbPyZ3H777R7r3NDMrHfv3h7r2Cm3h5vuX/Tiiy+GfrqnVL4m6/ldvHixx2eeeWbop6Xlm+tay0obAAAAAACAAuKhDQAAAAAAQAHVTHpUTsk477zzPM5L63R5sZb40qVv+Gy5ZOzhhx/usZYszWbOnOnxvffe63FRSsnm71J95eJKva4lmq6j6XD5PE6ZMqXJjqna6Hcxl5bV5Zv5e/jSSy95XOkl13mJuKbcnHjiiR5r2pRZPN5Zs2aFtlq69s6fPz+8njZtmsdf//rXQ9s3v/lNjw855BCPtSS3WUyV01jTAcxiakhOCdCyt7k8p2JZeGXk+8zxxx/vsZaXzjSFp9x5wubT73b+nn/00Uce6+eeU9T0vOq1smfPnqHf7rvv7rGOWbOYHqVpIvhset5ySffDDjvMY01tyulp/fr181hTbjQ9wyymbuS0V0350FSs/v37h36LFi3yuBpT8/Pc5KKLLvJ40qRJHg8ePDj022uvvTzOv+E01VBTFfP9Ts+PphTntGEdi1qC3czsL3/5S8njXbhwYein85tqOG8Nkc+1piLldGxNg9LPS6+zZvEc6nnTuZGZ2f777+9xvnbrGDvuuOM8fuWVV0K/Ipw3VtoAAAAAAAAUEA9tAAAAAAAACoiHNgAAAAAAAAVUM3vaaPkwM7Nx48Z5nPMXNUcx7y+A8jT3+thjjw1tP/zhDz3u0KGDxzkn++677/ZYy/A19f4wmj+spXbbtWsX+ul+Gx988EFoq4USjfU56aSTPNZcYi0xbUZO/pbQHNtczlL3Wso5vHo+NHc/l7jdVLo3Q5cuXULb9ddf77Huk5Lz83UfiN/97nehLecxV7N8DsaPH+/xihUrQpvmX3ft2tXjXMq2vvfXvY1yW97HQa97mkeer3m1vI9XJeWys2PHjvVY77N5bF922WUeUwK6ePSardfNIUOGhH56Hc17KTz44IMe63U/7ydWbizq9Vffvwj7NjQmHRN/+MMfQpveg3RPjDznW7Nmjce6B+arr74a+uk+c/m+qHNKLUud92d56KGHPK7G+2D+vumcWUurz5s3L/TTvfJ+9KMfhbYBAwZ4rOdU5zq5Tct/Dx8+PPTT85N/L/7gBz/wWMeU/o40i9+ZvC+OjtNaGov6/62fT0PpOP3FL34R2nRPvryv5ujRoz1esGCBx0X8/FlpAwAAAAAAUEA8tAEAAAAAACigqk6P2mqrrTy+7bbbQpsu/c5L1a666iqPKV+6eXQ5qJb3NTPr06ePx7rsTMtXmtW/bDSn0ejSurwsuL5+ecm5/ru87FHTDc4++2yPO3XqFPr99a9/9Tin0+mSzmpPG8if7emnn16y36233hpes4S/MnR5r1m8dmn5RDOzffbZx2MtFa3fZbOYBqNLf3PqzMiRIz2+4IIL6v1b+h55PMyZM8fjl19+ObQVcZlqY8mfy7JlyzzW1Bczs4cfftjjPfbYw+N169aFfnov1CXnb7zxRujXsWNHj0eNGhXazjjjDI/r6uo8Xr16dejHeK4MvfeZme2yyy4l++n3w2zjMqUoFr0G6nzplFNOCf00BW7lypWh7d577/W4XAp2TkFVOveppTGr9xK9jpmZ3XDDDR7rb4Y8N9Q5qn52Oe1Cz+9rr70W2k477TSPdWxrier8Ov9WqfbzVi79e+nSpR7reTOLWzNoSln+naDn/xvf+IbHuTR4fcdkFq/TmpaV58NaXj7PszTFuL5UKZSmc1vdWkN/v5nFz3LatGmhbfr06R4X/XcaK20AAAAAAAAKiIc2AAAAAAAABVTV6VFf/epXPc7VozR14Nxzzw1tOV0Hm06X9Obd8nW5oC7b7dGjR+inSxs7d+7s8dy5c0M/XZafl07qcehS1rxcuH///h4fddRRoW3vvfcu+R763mYx7SunR5188slWKzS1wixWCtI0G618gcrJlYA0nVDPhZlZmzZtPL7kkks8zt/tt956y+NDDz3U48MPPzz00/fPFRp03Oty7uXLl4d+WiWpGqtkNJQu183L7x955BGP85LfhtBqGs8991y9x6FVrHIqGypD093MYoqbLvW+7rrrQj9SuotN5yA6v9lzzz1DP71WTp06NbTp3Ef75XQK/Vt57kMaxsb/3zpP0XTjcmlmKvfT+2dObdL5Zfv27T3WlCqzOI9evHhxaKv29Khy9FzlNN877rjDY5275/Qo/Y3Yq1cvj/P8o1wlTn3Pfv36eZyrKi5cuLDkMZnFCrS1OhYbSueNmoqfaXWq888/P7S1pCq/rLQBAAAAAAAoIB7aAAAAAAAAFBAPbQAAAAAAAAqo6va00bxvLd2dy689//zzHk+ePDm0Fb3kV5HpnghTpkwJbZqbqyUUcznFgw8+2OPDDjvM43xeNLc05/ZqTqrmFee9b/SYdJ8Ps/hd0r+dSyuvX7/e4xtvvDG01VLOsZ43s1iKT78Xuh8GKkc/Y7NYHvqaa64Jbbr/UPfu3T3O+2PssMMOHut4yLnhKudka76wXnfHjRsX+s2aNcvjWho3W6LS+e96ndP9NsziPVSvr3mfHVRG3odP6X4O99xzT1McDipEx6zuh5fnQbrPhe6HYRbPv75fuesBe2U03KZ+drmfvn733XdD29q1a0u+Ry75rfsuzpgxY5OOo9bkz1bn5J06dfJYx5uZ2aBBgzzW3xNaTtws7jm00047hbZtttnGY50v5d+cOh/Ov2V0L8Fye1TBbPfddw+vzz77bI/1M897Co0ePdrjZcuWNdLRNT5W2gAAAAAAABQQD20AAAAAAAAKqOrSo4YPH+5x3759Pc4lvTQNgOXdlaNLDPOybU2f6dq1q8e6bNAspl5oCcW83LB169Yl/41ZXH6oSxZzKT9dkpzbdHmyLrd84YUXQr8JEyZ4/PTTT1st0c/9jDPOqLffzJkzPc4paqiMvOT2/vvv93iXXXYJbeedd57Huny4Q4cOoV+5NCily3jz9fSWW27xWFO2cplOUqKan15jBwwYENp03M6dO9dj0okbx9ChQ8Nr/ZyXL1/usab/ovg0FeKggw4q+d/NYgrUunXrGvS3SK9oXvr55/vb/PnzPdZ5s85rzeI9OZeK1t81tXyuy5Vu1+vmwIEDQz8dV88995zHPXr0CP00VTinMdbn9ddfD68ffvjhkn83HyM2tt1223l8++23h7a2bdt6rOPo7rvvDv3092hLHiustAEAAAAAACggHtoAAAAAAAAUEA9tAAAAAAAACqjF72mje6OYmU2cONFjzQ19/PHHQ7+cF4fK0FxSzdk1MzvnnHM87tatm8f5HG6//fYev/POOx7nfNQhQ4Z4rPvWmJnV1dV5vGrVKo/zPhqaW5pLUet7LlmyxON//OMfod/KlSs9rrX9WnSM9ezZM7Rpfqnu9cPeJU1D92i66aabQtvq1as9Hjt2rMd77bVX6Kf7TWkecN4jTMd6Lhuue+toGcaWnFdcrTRfP49n3ddLr6+cx8rRPYVyaVO9t+o9R6+zKD4t6azX27ynnpYxzvdM3f9G9/1jLBZXHqdvv/22xzoPzXs36p42udy03k917lnr3wP9/9d71QMPPBD66fjr3bu3x7169Qr9dIzlOf57773n8fTp0z2++uqrQz/d4ybvYaPvz/42n9LP5Mwzz/S4X79+9fbT32ljxowJ/fT+2ZKx0gYAAAAAAKCAeGgDAAAAAABQQC0yPUqXD1511VWhTZcP6nLuSy65JPTLy/tRGbq0T5f3msUUmWeeecZjXd5mFssMl1v6rf1yuUz9jpRbFqfLjvNx1LfsOC89reWlqNtuu63HOfVsm2228VjTE2v582ouuoTXLJZDfOqppzw+4IADQj9dJrxgwQKPFy1aFPq98sorHn/44YehjeW+LYdeU3Mp6Z133tljTeXIy/mrZRlyc9CSvrm0rM5ZdBk446vYdEyZxTRvTcHO40bHVfv27UObznf0fprnMOVwH25eb775psc6Hx45cmS9/2bAgAHh9bJlyyp/YFVGr5svvvhiaNPr6Lx58zx+9hdvJsMAAANESURBVNlnQ7/u3bt7nMfzzJkzPV64cKHHuax3ufHGWNyYpmePGzfO4/xbT+e23//+9z3W9MNqwkobAAAAAACAAuKhDQAAAAAAQAG1yPSoTp06eTxs2LDQpukuWtFEU3PQ/MqlG23qcm/tV4mqTSxR3Hxr1671+KCDDgpturybCifFUl+qhcaoPZr29M9//jO06bLw5557zmPScyrn/fff91irupmZHXXUUR7ruaEaX7HleYUu51+zZo3H+R6pc5o8f9XUc+YtLZN+D2bNmuVxHs+aepy3HODau3ny56XXW/2cNUbT2GqrrcLra665xuM2bdp4nK+Td9xxh8ezZ89upKMrDlbaAAAAAAAAFBAPbQAAAAAAAAqIhzYAAAAAAAAF1Gpz8mFbtWrVbMmzWsrw3HPP9Xj8+PGhn+bF/epXv/L4iiuuCP1aWi7ohg0bNr2WYxnNeQ5hszds2LBfJd6I89h8GItVgbH4GXLJaS1HreU0m/NeylisCjU7FnW+mksJ69w8lwMv4j42jMXNo+e7devWHuuenWZx75u8d6PucVOhva1qdixWk5Y4FrXEt5nZk08+6XHbtm091r1qzcxGjhzp8apVqxrp6JpFybHIShsAAAAAAIAC4qENAAAAAABAAbWYkt+6HHTx4sUer1u3LvSbMWOGx1oKrKWlQwEA0FzWr19f9jWALZPTnlA79DeJpkBpGWqz+Nsnp9ABLZmmhx5zzDH1tq1Zs8bjCRMmhH5r165tpKMrJq4AAAAAAAAABcRDGwAAAAAAgALioQ0AAAAAAEABtZiS3/XRUuBZEcsiNlRLLOGGjVBOsQowFqsCY7EKMBarAmOxCjAWqwJjsQpUw1jceuutPc6l7msEJb8BAAAAAABaCh7aAAAAAAAAFNDmlvyuM7OljXEgDVVNKVBl7FbB9yrcOawhnMeWj3NYHTiPLR/nsDpwHls+zmF14Dy2fFVxDms0JUqVPI+btacNAAAAAAAAmgbpUQAAAAAAAAXEQxsAAAAAAIAC4qENAAAAAABAAfHQBgAAAAAAoIB4aAMAAAAAAFBAPLQBAAAAAAAoIB7aAAAAAAAAFBAPbQAAAAAAAAqIhzYAAAAAAAAF9H+940acpHi3zQAAAABJRU5ErkJggg==\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: 3.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(4*1,4)) # encoded 64 # (8,8) # (4, 4*4)\n",
    "    #plt.imshow(decoded_imgs[i].reshape(4, 14*14).T) # input or decoded 784 # (28,28) # (4, 14*14)\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "n = 10\n",
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(decoded_imgs[i].reshape(28,28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE = 0.0042, PSNR (dB) = 23.77\n",
      "MSE = 0.0157, PSNR (dB) = 18.04\n",
      "MSE = 0.0014, PSNR (dB) = 28.62\n",
      "MSE = 0.0123, PSNR (dB) = 19.10\n",
      "MSE = 0.0079, PSNR (dB) = 21.03\n",
      "MSE = 0.0007, PSNR (dB) = 31.55\n",
      "MSE = 0.0199, PSNR (dB) = 17.02\n",
      "MSE = 0.0247, PSNR (dB) = 16.06\n",
      "MSE = 0.0266, PSNR (dB) = 15.74\n",
      "MSE = 0.0126, PSNR (dB) = 19.01\n",
      "Wall time: 4.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": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9361\n",
      "0.8336\n",
      "0.9779\n",
      "0.8918\n",
      "0.8953\n",
      "0.9899\n",
      "0.7283\n",
      "0.6917\n",
      "0.7225\n",
      "0.8543\n",
      "Wall time: 399 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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9816\n",
      "0.9431\n",
      "0.9947\n",
      "0.9555\n",
      "0.9575\n",
      "0.9967\n",
      "0.9195\n",
      "0.8979\n",
      "0.8660\n",
      "0.9565\n",
      "Wall time: 890 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": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.save(\"./Deep Autoencoders in Keras - 784-256-128-64-32-16-32-64-128-256-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
}