vdsr4dem_upscale_3.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "name": "vdsr4dem_upscale_3.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "machine_shape": "hm"
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "idupnGQ-I-UC"
      },
      "source": [
        "VDSR for DEM processing - Upscale factor 3\n",
        "====================================\n",
        "Very Deep Super Resolution (VDSR) CNN based on the implementation of George Seif [GitHub](https://github.com/GeorgeSeif/VDSR-Keras).\n",
        "Thanks to [Zushicat](https://github.com/zushicat) for the support.\n",
        "\n",
        "An implementation of the Super Resolution CNN proposed in:\n",
        "\n",
        "[Kim, Jiwon, Jung Kwon Lee, and Kyoung Mu Lee. \"Accurate image super-resolution using very deep convolutional networks.\" Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. 2016.](https://arxiv.org/abs/1511.04587)\n",
        "\n",
        "This Notebook is running on Google Colab with data storage on Google Drive. Of course this could be adapted to local resources.\n",
        "\n",
        "For background information and results visit: [https://opendem.info/superResolution.html](https://opendem.info/superResolution.html).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "xXvhPyws9-T7",
        "outputId": "8e5e0d3f-f21b-4022-cfe1-0113d1f75bd4"
      },
      "source": [
        "import os\n",
        "# connect with google drive  \n",
        "from google.colab import drive\n",
        "drive.mount(\"/content/drive\")"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Mounted at /content/drive\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gIqPFMfvI9AE"
      },
      "source": [
        "Training\n",
        "--------\n",
        "**This part could be skipped if you like to use the prepared weights.**\n",
        "\n",
        "Download the training, test and validation data of the austrian DTM [here](https://www.openmaps.online/data/austria_dtm.zip). License: Data by - Geoland.at (2019) - CC BY 4.0 (CC BY 4.0)\n",
        "\n",
        "When using own data, make sure that the low resolution images are upsampled to 300*300 px. Input size for both datasets must be 300 * 300 px for this CNN."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3_AGZx2bNS8g"
      },
      "source": [
        "Get rasterio for GeoTiff processing"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vxoYac2yx9lS",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "0a74c57a-39bd-4f22-b1ff-0c6aebd07c61"
      },
      "source": [
        "!pip install rasterio"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pPMTTGgtZHYd",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "4cc326e1-77ae-452e-e387-27c8385d6f77"
      },
      "source": [
        "import os\n",
        "import tensorflow as tf\n",
        "import numpy as np\n",
        "import rasterio\n",
        "from tensorflow.python.data.experimental import AUTOTUNE\n",
        "from rasterio.plot import reshape_as_image, reshape_as_raster, show\n",
        "import re\n",
        "\n",
        "# adapt your paths\n",
        "HR_PATH_TRAIN = \"/content/drive/My Drive/vdsr4dem/hr_train\"\n",
        "LR_PATH_TRAIN = \"/content/drive/My Drive/vdsr4dem/lr_train\"\n",
        "HR_PATH_VAL = \"/content/drive/My Drive/vdsr4dem/hr_valid\"\n",
        "LR_PATH_VAL = \"/content/drive/My Drive/vdsr4dem/lr_valid\"\n",
        "\n",
        "# maximal and minimal values of the whole trainings dataset\n",
        "ds_max = 3150\n",
        "ds_min = 170\n",
        "\n",
        "# use human sorting to avoid problems with the order of the files for training and validation\n",
        "def atoi(text):\n",
        "    return int(text) if text.isdigit() else text\n",
        "\n",
        "def natural_keys(text):\n",
        "    return [ atoi(c) for c in re.split(r'(\\d+)', text) ]\n",
        "\n",
        "def _images_dataset(image_files):    \n",
        "    listAll = []\n",
        "    for img in image_files:\n",
        "        with rasterio.open(img) as src:\n",
        "            array = src.read(\n",
        "                out_shape=(\n",
        "                    src.count,\n",
        "                    int(src.height),\n",
        "                    int(src.width)\n",
        "                )\n",
        "            )\n",
        "            array = reshape_as_image(array)            \n",
        "            lr = (array - ds_min) / (ds_max - ds_min)\n",
        "            listAll.append(lr)        \n",
        "    ds = tf.data.Dataset.from_tensor_slices(listAll)\n",
        "    #print(ds)\n",
        "    return ds\n",
        "\n",
        "def _get_dataset(path, num_images):\n",
        "    file_names = [f\"{path}/{x}\" for x in os.listdir(path) if x.find(\".tif\") != -1]\n",
        "    if num_images is not None:\n",
        "        file_names = file_names[:num_images]\n",
        "    print(\"--->\", len(file_names))\n",
        "    file_names.sort(key=natural_keys)\n",
        "    print(file_names)\n",
        "    ds = _images_dataset(file_names)    \n",
        "    return ds\n",
        "\n",
        "\n",
        "def dataset(batch_size=16, repeat_count=None, random_transform=True, use_set=\"train\", num_images=None):\n",
        "    print(num_images)\n",
        "    if use_set == \"train\":\n",
        "        print('LR_PATH_TRAIN')\n",
        "        lr_ds = _get_dataset(LR_PATH_TRAIN, num_images)\n",
        "        print('HR_PATH_TRAIN')\n",
        "        hr_ds = _get_dataset(HR_PATH_TRAIN, num_images)\n",
        "    elif use_set == \"val\":\n",
        "        print('LR_PATH_VAL')\n",
        "        lr_ds = _get_dataset(LR_PATH_VAL, num_images)\n",
        "        print('HR_PATH_VAL')\n",
        "        hr_ds = _get_dataset(HR_PATH_VAL, num_images)\n",
        "    else:\n",
        "        return None\n",
        "    \n",
        "    ds = tf.data.Dataset.zip((lr_ds, hr_ds))\n",
        "    ds = ds.batch(batch_size)\n",
        "    ds = ds.prefetch(buffer_size=AUTOTUNE)\n",
        "    return ds\n",
        "\n",
        "if __name__ == \"__main__\":\n",
        "    # To avoid problems with the order of the files always set num_images to the real number of images in the folder!\n",
        "    train_ds = dataset(batch_size=16, use_set=\"train\", num_images=30)\n",
        "    val_ds = dataset(batch_size=1, use_set=\"val\", num_images=30)\n",
        "    print(\"done all\")"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6twCsQBWJkXY"
      },
      "source": [
        "VDSR Model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jzoWQQHCnksk"
      },
      "source": [
        "from keras import Sequential\n",
        "from keras.layers import Conv2D\n",
        "from keras.models import Sequential, Model\n",
        "from keras.layers import Dense, Activation, Input, add\n",
        "from keras import optimizers\n",
        "from keras import losses\n",
        "from keras.optimizers import Adam\n",
        "import tensorflow as tf\n",
        "\n",
        "IMG_SIZE = (300, 300, 1)\n",
        "\n",
        "model = Sequential()\n",
        "\n",
        "input_img = Input(shape=IMG_SIZE)\n",
        "\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(input_img)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "\n",
        "model = Conv2D(1, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "\n",
        "res_img = model\n",
        "\n",
        "output_img = add([res_img, input_img])\n",
        "\n",
        "model = Model(input_img, output_img)\n",
        "\n",
        "adam = Adam(lr=0.0001)\n",
        "model.compile(adam, loss='mse', metrics=[tf.keras.metrics.MeanAbsoluteError()])"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8Ue_jX8wJ46A"
      },
      "source": [
        "Fit the Model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "34jeS3_Y8uop",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "6930f2a7-b273-454d-e68d-786f1c92b33f"
      },
      "source": [
        "from keras.callbacks import ModelCheckpoint\n",
        "\n",
        "# adapt path\n",
        "checkpoint = ModelCheckpoint(\"/content/drive/My Drive/vdsr4dem/weights/vdsr4dem_upscale_3.h5\", monitor='loss', verbose=1, save_best_only=True, mode='auto', save_weights_only=True)\n",
        "# adapt epochs\n",
        "model.fit(train_ds, epochs=25, validation_data=val_ds, callbacks=[checkpoint])"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<tensorflow.python.keras.callbacks.History at 0x7f1b98242a20>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xw3dj85BNTRM"
      },
      "source": [
        "Visualization of the loss"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        },
        "id": "NzDiuhbEFIAX",
        "outputId": "fccec8f3-dee2-4398-dfa3-347f7c304bb3"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "plt.plot(model.history.history['loss'])\n",
        "plt.title('model loss')\n",
        "plt.ylabel('loss')\n",
        "plt.xlabel('epoch')\n",
        "plt.legend(['train', 'test'], loc='upper left')\n",
        "plt.show()"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df5RdZX3v8fdnfs8kk5khDAgJmFSwJmANZhKxFK5XLzSiEiq/gqggLNBVWeq1dYltxV6u3ou37dVrS5VQIj8KBAqiUUNTqeJtbwUzhAgJSAlpJBPQhJAJgfycme/9Y++ZnJycyZwzZ3bOZM7ntdZZs8+zn73PszMr+WQ/zz7Po4jAzMxstGoq3QAzMzuyOUjMzKwsDhIzMyuLg8TMzMriIDEzs7I4SMzMrCwOErMMSbpN0peLrLtB0n8p9zxmh5uDxMzMyuIgMTOzsjhIrOqlXUqfk/SkpNcl3SrpWEkPSdoh6WFJHTn1z5O0VlKvpEckzcrZd5qkVelx9wJNeZ/1fkmr02P/TdLvjLLNV0taJ+kVScskHZ+WS9LXJG2W9KqkpySdmu47V9LTads2SfrjUf2BmeVxkJglLgDOBt4MfAB4CPgToJPk78mnACS9GbgH+Ey6bznwfUkNkhqA7wJ3AkcB/5Cel/TY04AlwMeBqcDNwDJJjaU0VNK7gf8JXAwcB/wKWJruPgc4K72OtrTO1nTfrcDHI6IVOBX4cSmfazYcB4lZ4q8j4jcRsQn4F+CxiHgiInYDDwKnpfUuAX4YET+KiH3AXwLNwO8CpwP1wNcjYl9E3A+szPmMa4CbI+KxiOiPiNuBPelxpbgMWBIRqyJiD/AF4J2SZgD7gFbgLYAi4pmIeCk9bh8wW9KUiNgWEatK/FyzghwkZonf5GzvKvB+crp9PMkdAAARMQBsBKal+zbFgTOh/ipn+43AH6XdWr2SeoET0uNKkd+G10juOqZFxI+BvwFuAjZLWixpSlr1AuBc4FeSfirpnSV+rllBDhKz0rxIEghAMiZBEgabgJeAaWnZoBNztjcCX4mI9pxXS0TcU2YbJpF0lW0CiIhvRMRcYDZJF9fn0vKVEbEQOIakC+6+Ej/XrCAHiVlp7gPeJ+k9kuqBPyLpnvo34GdAH/ApSfWSPgjMzzn2FuATkt6RDopPkvQ+Sa0ltuEe4GOS5qTjK/+DpCtug6R56fnrgdeB3cBAOoZzmaS2tEvuVWCgjD8HsyEOErMSRMSzwIeBvwZeJhmY/0BE7I2IvcAHgSuAV0jGU76Tc2w3cDVJ19M2YF1at9Q2PAx8EXiA5C7oTcCidPcUksDaRtL9tRX4i3TfR4ANkl4FPkEy1mJWNnlhKzMzK4fvSMzMrCwOEjMzK4uDxMzMyuIgMTOzstRVugGHw9FHHx0zZsyodDPMzI4ojz/++MsR0TlSvaoIkhkzZtDd3V3pZpiZHVEk/WrkWu7aMjOzMjlIzMysLA4SMzMrS1WMkRSyb98+enp62L17d6WbkqmmpiamT59OfX19pZtiZhNU1QZJT08Pra2tzJgxgwMna504IoKtW7fS09PDzJkzK90cM5ugMu3akrRA0rPpkqDXFdh/VrosaZ+kC3PK/3O6HOnga7ek89N9t0n6j5x9c0bTtt27dzN16tQJGyIAkpg6deqEv+sys8rK7I5EUi3J4jpnAz3ASknLIuLpnGovkMx+esDa0RHxE2BOep6jSGZJ/aecKp9LV58rt43lnmLcq4ZrNLPKyrJraz6wLiLWA0haCiwEhoIkIjak+w61LsKFwEMRsTO7pha2bedeBgaCqZNLWlLbzKyqZNm1NY1kRbhBPWlZqRaRLOST6yuSnpT0tXRhn4NIukZSt6TuLVu2jOJjYfvOfWx9fe+ojh1Jb28vf/u3f1vyceeeey69vb0ZtMjMbHTG9eO/ko4D3gqsyCn+AvAWYB5wFPD5QsdGxOKI6IqIrs7OEb/hX1BtjegfyGa9luGCpK+v75DHLV++nPb29kzaZGY2GlkGySaStawHTU/LSnEx8GC6NCgAEfFSJPYA3+bApUzHVJZBct111/H8888zZ84c5s2bx5lnnsl5553H7NmzATj//POZO3cup5xyCosXLx46bsaMGbz88sts2LCBWbNmcfXVV3PKKadwzjnnsGvXrkzaamZ2KFmOkawETpY0kyRAFgEfKvEcl5LcgQyRdFxEvKRkFPl8YE25Df1v31/L0y++elD5vv4B9vYNMKmx9D+m2cdP4UsfOGXY/TfeeCNr1qxh9erVPPLII7zvfe9jzZo1Q4/pLlmyhKOOOopdu3Yxb948LrjgAqZOnXrAOZ577jnuuecebrnlFi6++GIeeOABPvzhD5fcVjOzcmR2RxIRfcC1JN1SzwD3RcRaSTdIOg9A0jxJPcBFwM2S1g4eL2kGyR3NT/NOfZekp4CngKOBL2d1DYPPOx2OxYjnz59/wHc9vvGNb/C2t72N008/nY0bN/Lcc88ddMzMmTOZMyd5+nnu3Lls2LDhMLTUzOxAmX4hMSKWA8vzyq7P2V5J0uVV6NgNFBicj4h3j20rGfbOoXfnXl54ZSdvPraVpvrasf7YA0yaNGlo+5FHHuHhhx/mZz/7GS0tLbzrXe8q+F2Qxsb9zxnU1ta6a8vMKmJcD7ZXWm1Nck/Sl8E4SWtrKzt27Ci4b/v27XR0dNDS0sIvf/lLHn300TH/fDOzsVK1U6QUoy4NkiwG3KdOncoZZ5zBqaeeSnNzM8cee+zQvgULFvCtb32LWbNm8du//ducfvrpY/75ZmZjRRGHYwSgsrq6uiJ/YatnnnmGWbNmHfK4vX39/PLXO5je0cJRkxqybGKmirlWM7N8kh6PiK6R6rlr6xBqa5I/nv6BQ33x3sysujlIDqFGIJTJGImZ2URR1UEyUreepEy/lHg4VEPXpZlVVtUGSVNTE1u3bh3xH9ojOUgG1yNpamqqdFPMbAKr2qe2pk+fTk9PDyNN6Lhlxx4E7Nx8ZM4APLhCoplZVqo2SOrr64taNfCvbl/Jpt7dPPTpUa2fZWY24VVt11ax2lsa2L4zm6nkzcwmAgfJCNqb6+ndtW/kimZmVcpBMoL2lnp27u1nT19/pZtiZjYuOUhG0NaSfKN9+07flZiZFeIgGUFHSz2Au7fMzIbhIBlBe3NyR9LrOxIzs4IcJCNoH7wj8ZNbZmYFOUhG0NY8GCS+IzEzK8RBMoKOdPr43l2+IzEzK8RBMoJJDbXU1ch3JGZmw8g0SCQtkPSspHWSriuw/yxJqyT1Sbowb1+/pNXpa1lO+UxJj6XnvFdSpitOSaK9xV9KNDMbTmZBIqkWuAl4LzAbuFTS7LxqLwBXAHcXOMWuiJiTvs7LKf8q8LWIOAnYBlw15o3P09Zc78F2M7NhZHlHMh9YFxHrI2IvsBRYmFshIjZExJNAUUsQShLwbuD+tOh24Pyxa3JhHS0N7toyMxtGlkEyDdiY874nLStWk6RuSY9KGgyLqUBvRPSNdE5J16THd480VfxI2lvqHSRmZsMYz4Ptb0wXnf8Q8HVJbyrl4IhYHBFdEdHV2dlZVkPamhvY7jESM7OCsgySTcAJOe+np2VFiYhN6c/1wCPAacBWoF3S4DoqJZ1ztNpb6tnmMRIzs4KyDJKVwMnpU1YNwCJg2QjHACCpQ1Jjun00cAbwdCTr4v4EGHzC63Lge2Pe8jwdngHYzGxYmQVJOo5xLbACeAa4LyLWSrpB0nkAkuZJ6gEuAm6WtDY9fBbQLekXJMFxY0Q8ne77PPBZSetIxkxuzeoaBg3NAOzuLTOzg2S61G5ELAeW55Vdn7O9kqR7Kv+4fwPeOsw515M8EXbYtKfTpGzfuY9jWpsO50ebmY1743mwfdzoSO9ItvnJLTOzgzhIiuAZgM3MhucgKcLQDMAeIzEzO4iDpAiDdyRebtfM7GAOkiJMbqyjrkb+LomZWQEOkiJ4BmAzs+E5SIrU1lzvri0zswIcJEVqb2lw15aZWQEOkiJ1eAZgM7OCHCRF8gzAZmaFOUiKlKxJ4q4tM7N8DpIitTfX8/refvb2FbWYo5lZ1XCQFKl9UjLfVu8u35WYmeVykBQpdwZgMzPbz0FSpKGJGz3gbmZ2AAdJkdqb06nkX3fXlplZLgdJkXxHYmZWmIOkSJ4B2MysMAdJkSY31lFbIz+1ZWaWx0FSJEm0N9d7uV0zszyZBomkBZKelbRO0nUF9p8laZWkPkkX5pTPkfQzSWslPSnpkpx9t0n6D0mr09ecLK8hV3uLZwA2M8tXl9WJJdUCNwFnAz3ASknLIuLpnGovAFcAf5x3+E7goxHxnKTjgcclrYiI3nT/5yLi/qzaPpz2lgZ3bZmZ5cksSID5wLqIWA8gaSmwEBgKkojYkO47YN6RiPj3nO0XJW0GOoFeKqi9uZ5fv7q7kk0wMxt3suzamgZszHnfk5aVRNJ8oAF4Pqf4K2mX19ckNQ5z3DWSuiV1b9mypdSPLajNU8mbmR1kXA+2SzoOuBP4WEQM3rV8AXgLMA84Cvh8oWMjYnFEdEVEV2dn55i0p6OlwTMAm5nlyTJINgEn5LyfnpYVRdIU4IfAn0bEo4PlEfFSJPYA3ybpQjssPAOwmdnBsgySlcDJkmZKagAWAcuKOTCt/yBwR/6genqXgiQB5wNrxrTVhzD0pUR/u93MbEhmQRIRfcC1wArgGeC+iFgr6QZJ5wFImiepB7gIuFnS2vTwi4GzgCsKPOZ7l6SngKeAo4EvZ3UN+dpa0qnk3b1lZjYky6e2iIjlwPK8sutztleSdHnlH/f3wN8Pc853j3Ezi9bh+bbMzA4yrgfbx5vBGYD95JaZ2X4OkhIMzQDsri0zsyEOkhLsDxLfkZiZDXKQlMAzAJuZHcxBUoLBGYB9R2Jmtp+DpERtLfV+asvMLIeDpESeJsXM7EAOkhK5a8vM7EAOkhJ5BmAzswM5SErU3uyuLTOzXA6SEnW0eAZgM7NcDpISeQZgM7MDOUhKNDgD8HZ/KdHMDHCQlKy9Obkj2eYBdzMzwEFSso4WzwBsZpbLQVIizwBsZnYgB0mJ2jzYbmZ2AAdJiVrTGYC3+Y7EzAxwkJTMMwCbmR0o0yCRtEDSs5LWSbquwP6zJK2S1Cfpwrx9l0t6Ln1dnlM+V9JT6Tm/IUlZXkMhngHYzGy/zIJEUi1wE/BeYDZwqaTZedVeAK4A7s479ijgS8A7gPnAlyR1pLu/CVwNnJy+FmR0CcNqb65nu+9IzMyAbO9I5gPrImJ9ROwFlgILcytExIaIeBLIn2/k94EfRcQrEbEN+BGwQNJxwJSIeDQiArgDOD/DayiovaXBYyRmZqksg2QasDHnfU9aVs6x09LtEc8p6RpJ3ZK6t2zZUnSji9HuGYDNzIZM2MH2iFgcEV0R0dXZ2Tmm525vbvDjv2ZmqSyDZBNwQs776WlZOcduSrdHc84x095Sz2t7+tjX7xmAzcyyDJKVwMmSZkpqABYBy4o8dgVwjqSOdJD9HGBFRLwEvCrp9PRprY8C38ui8Yey/9vtvisxM8ssSCKiD7iWJBSeAe6LiLWSbpB0HoCkeZJ6gIuAmyWtTY99BfjvJGG0ErghLQP4Q+DvgHXA88BDWV3DcNo9A7CZ2ZC6LE8eEcuB5Xll1+dsr+TArqrcekuAJQXKu4FTx7alpRmcAdh3JGZmE3iwPUvu2jIz289BMgqDU8n7uyRmZkUGiaRPS5qixK3ptCbnZN248cozAJuZ7VfsHcmVEfEqydNTHcBHgBsza9U4NzgDsLu2zMyKD5LBiRHPBe6MiLU5ZVVHEm3N9fT6qS0zs6KD5HFJ/0QSJCsktXLw/FhVpb2l3uu2m5lR/OO/VwFzgPURsTOdnfdj2TVr/PMMwGZmiWLvSN4JPBsRvZI+DPwZsD27Zo1/7S0N7toyM6P4IPkmsFPS24A/IvlG+R2ZteoI4FUSzcwSxQZJX7r+x0LgbyLiJqA1u2aNf+0tDQ4SMzOKHyPZIekLJI/9nimpBqjPrlnjX+4MwPW1/l6nmVWvYv8FvATYQ/J9kl+TzI/1F5m16gjQ7i8lmpkBRQZJGh53AW2S3g/sjoiqHiNp88SNZmZA8VOkXAz8nGS694uBxyRdmGXDxrvB+bZ6Pd+WmVW5YsdI/hSYFxGbASR1Ag8D92fVsPHOMwCbmSWKHSOpGQyR1NYSjp2Q2pvTOxKPkZhZlSv2juQfJa0A7knfX0LeglXVpm3ojsRdW2ZW3YoKkoj4nKQLgDPSosUR8WB2zRr/pjR5BmAzMyhhqd2IeAB4IMO2HFE8A7CZWeKQ4xySdkh6tcBrh6RXRzq5pAWSnpW0TtJ1BfY3Sro33f+YpBlp+WWSVue8BiTNSfc9kp5zcN8xo7v08nmaFDOzEe5IImLU06BIqgVuAs4GeoCVkpZFxNM51a4CtkXESZIWAV8FLomIu0i+t4KktwLfjYjVOcddFhHdo23bWGlrcZCYmWX55NV8YF1ErI+IvcBSkrm6ci0Ebk+37wfeIyl/waxL02PHnQ7PAGxmlmmQTAM25rzvScsK1omIPpKp6afm1bmE/U+LDfp22q31xQLBA4CkayR1S+resmXLaK/hkNy1ZWY2zr8LIukdwM6IWJNTfFlEvBU4M319pNCxEbE4IroioquzszOT9rW1eHErM7Msg2QTcELO++lpWcE6kuqANpIvOw5aRN7dSERsSn/uAO4m6UKriPbmBnakMwCbmVWrLINkJXCypJmSGkhCYVlenWXA5en2hcCP03VPSKeqv5ic8RFJdZKOTrfrgfcDa6iQjkmeAdjMrOjvkZQqIvokXQusAGqBJRGxVtINQHdELANuBe6UtA54hSRsBp0FbIyI9TlljcCKNERqSeb7uiWraxhJ7gzAR09urFQzzMwqKrMgAYiI5eRNpRIR1+ds7yaZUbjQsY8Ap+eVvQ7MHfOGjlJ7OgPwdj+5ZWZVbFwPto93Hel8W9ted9eWmVUvB0kZPAOwmZmDpCyeAdjMzEFSltbGOmrkp7bMrLo5SMpQUyPaWxrY5jsSM6tiDpIyeZoUM6t2DpIytbXUu2vLzKqag6RMviMxs2rnIClTh8dIzKzKOUjK5BmAzazaOUjK5BmAzazaOUjK1J5+KfFVD7ibWZVykJRpMEi2uXvLzKqUg6RMngHYzKqdg6RM7TlrkpiZVSMHSZnaWxwkZlbdHCRlGuza8ndJzKxaOUjK5BmAzazaOUjKVFMj2jxNiplVsUyDRNICSc9KWifpugL7GyXdm+5/TNKMtHyGpF2SVqevb+UcM1fSU+kx35CkLK+hGJ5K3syqWWZBIqkWuAl4LzAbuFTS7LxqVwHbIuIk4GvAV3P2PR8Rc9LXJ3LKvwlcDZycvhZkdQ3FavcMwGZWxbK8I5kPrIuI9RGxF1gKLMyrsxC4Pd2+H3jPoe4wJB0HTImIRyMigDuA88e+6aXxDMBmVs2yDJJpwMac9z1pWcE6EdEHbAempvtmSnpC0k8lnZlTv2eEcwIg6RpJ3ZK6t2zZUt6VjKC9pYFefyHRzKrUeB1sfwk4MSJOAz4L3C1pSikniIjFEdEVEV2dnZ2ZNHJQW3M9va/7jsTMqlOWQbIJOCHn/fS0rGAdSXVAG7A1IvZExFaAiHgceB54c1p/+gjnPOw6WjwDsJlVryyDZCVwsqSZkhqARcCyvDrLgMvT7QuBH0dESOpMB+uR9Fskg+rrI+Il4FVJp6djKR8FvpfhNRTFMwCbWTWry+rEEdEn6VpgBVALLImItZJuALojYhlwK3CnpHXAKyRhA3AWcIOkfcAA8ImIeCXd94fAbUAz8FD6qqihaVJ27WPq5MYKt8bM7PDKLEgAImI5sDyv7Pqc7d3ARQWOewB4YJhzdgOnjm1LyzM4TUqvv0tiZlVovA62H1E8A7CZVTMHyRjwDMBmVs0cJGOgvTnt2vJgu5lVIQfJGGhtSmYA9hiJmVUjB8kY8AzAZlbNHCRjJJkmxUFiZtXHQTJGkjsSd22ZWfVxkIyRjhZ3bZlZdXKQjBHPAGxm1cpBMkY82G5m1cpBMkbaW+rZsbuPPs8AbGZVxkEyRjrS+ba85K6ZVRsHyRjJnQHYzKyaOEjGSJsnbjSzKuUgGSPtQ11bfnLLzKqLg2SMdKRdW9u8druZVRkHyRjxDMBmVq0cJGOktakOCbZ7mhQzqzIOkjEyNAOw70jMrMo4SMZQR0sD2/zUlplVmUyDRNICSc9KWifpugL7GyXdm+5/TNKMtPxsSY9Leir9+e6cYx5Jz7k6fR2T5TWUwjMAm1k1qsvqxJJqgZuAs4EeYKWkZRHxdE61q4BtEXGSpEXAV4FLgJeBD0TEi5JOBVYA03KOuywiurNq+2i1t9TzyusOEjOrLlnekcwH1kXE+ojYCywFFubVWQjcnm7fD7xHkiLiiYh4MS1fCzRLasywrWOi3RM3mlkVyjJIpgEbc973cOBdxQF1IqIP2A5MzatzAbAqIvbklH077db6oiQV+nBJ10jqltS9ZcuWcq6jaO0tDWxz15aZVZlxPdgu6RSS7q6P5xRfFhFvBc5MXx8pdGxELI6Irojo6uzszL6xeAZgM6tOWQbJJuCEnPfT07KCdSTVAW3A1vT9dOBB4KMR8fzgARGxKf25A7ibpAttXGhP59t6dXdfhVtiZnb4ZBkkK4GTJc2U1AAsApbl1VkGXJ5uXwj8OCJCUjvwQ+C6iPh/g5Ul1Uk6Ot2uB94PrMnwGkoyON+Wu7fMrJpkFiTpmMe1JE9cPQPcFxFrJd0g6by02q3AVEnrgM8Cg48IXwucBFyf95hvI7BC0pPAapI7mluyuoZSDU0l7wF3M6simT3+CxARy4HleWXX52zvBi4qcNyXgS8Pc9q5Y9nGseQZgM2sGo3rwfYjTbvXJDGzKuQgGUODy+3+5NktbOrdVeHWmJkdHpl2bVWbKc11vPstx/D9X7zID558kdNnTuUP3j6N9576Blqb6ivdPDOzTCgiKt2GzHV1dUV39+GbUWXjKzt58IlNfGdVDxu27qSpvobfP+UNfPDt0znjTVOpq/WNoJmNf5Iej4iuEes5SLITEax6oZcHn+jh+794ie279tHZ2sj5c47ng2+fzqzjphz2NpmZFctBkqNSQZJrT18/P/nlZr6zahM/eXYz+/qDWcdN4YOnTWPhnOM5ZkpTRdtnZpbPQZJjPARJrlde38sPnnyRB1Zt4hcbe2moreGvP3Qav3/KGyrdNDOzIcUGiTvrK+CoSQ189J0z+N4nz+Dhz/4nZh8/hU/etYp/XPPrSjfNzKxkDpIKO+mYydxx1XzeOr2Na+9exT+ueanSTTIzK4mDZByY0lTPHVfO53emt3Ht3U/w0FMOEzM7cjhIxonWpnpuv3I+bzuhnWvvcZiY2ZHDQTKODIbJnDRMfvikw8TMxj8HyTgzubGO26+cz2kntPOppU/wgydfHPkgM7MKcpCMQ5Mb67jtyvm8/cR2Pr10Nd//hcPEzMYvB8k4Nbmxjm9/bD5zT+zgM/c6TMxs/HKQjGNJmMxj7hs7+PTSJ1jmMDGzcchBMs5Naqzj21fMo2vGUXxm6RN8b3X+svdmZpXlIDkCTGqs47aPzWPejKP4r/eu5rtPOEzMbPzweiRHiJaGpJvryttW8tn7VrNh6+u85Q1T6Gxt5JjWRjpbG2mqr610M82sCmUaJJIWAP8HqAX+LiJuzNvfCNxBsg77VuCSiNiQ7vsCcBXQD3wqIlYUc86JrKWhjiVXzOPjdz7O1x9+7qD9U5rqOGZKE52TGzlmSuPQz2Nam+hsbaS1qY7m+lqa6mtpbqiluT551dSoAldjZhNFZkEiqRa4CTgb6AFWSloWEU/nVLsK2BYRJ0laBHwVuETSbGARcApwPPCwpDenx4x0zgmtpaGOO66cz5Yde9i8Yw9b0tfmHbuH3m/esYcnXuhl847d7N43MOI5G+tqDgiWwaBpqq+hobaGhroaGupqh7Yb69KyoX0527U11NWKutoa6mpEXY2oHyyrqaE+Z199bQ21aZ0aCQlqakSN2P9eSl+gnJ9ifx1gqK7Yv19pXTPLVpZ3JPOBdRGxHkDSUmAhkPuP/kLgz9Pt+4G/UfI3fyGwNCL2AP8haV16Poo454QniWOmNI24hklE8NqevqGAeX1PH7v29bNzbz+79/Wza28/u/Ylr91D2wPsGty/r59Xd/Wxt2+Avf0D7O0bYE/fAHv7+ofeDxwBqxBIHBwwJIU6oI4OqpvszK23vw559ZLM2h9cQyGX15ak7OB6OR+X1/6Rw/Cgcxx0zsLnONSph9s1XHtGFdmHIedH8xET6T8gSy6fx4lTWzL9jCyDZBqwMed9D/CO4epERJ+k7cDUtPzRvGOnpdsjnRMASdcA1wCceOKJo7uCI5wkWpvqaW2q502dkzP5jL7+/SEzGDh9/UHfwAB9A0Fff7CvP9nel7NvX38MbQ9EMDAAAxFEJD8HYvD9/u2BIH2f1Av2HwMM1U32pXUiCPLKct4nBw6WHbg/2bX//IOfEUPb+88xeBw5+4ZOnld2wPk4MInzlwfKz+lCywfln2OEtznnGv5/AcMfU1r9QzkcayGN6hOOgP8claKhLvtnqibsYHtELAYWQ7KwVYWbM2HV1dZQV1tDS0OlW2JmlZJlVG0CTsh5Pz0tK1hHUh3QRjLoPtyxxZzTzMwOoyyDZCVwsqSZkhpIBs+X5dVZBlyebl8I/DiS+91lwCJJjZJmAicDPy/ynGZmdhhl1rWVjnlcC6wgeVR3SUSslXQD0B0Ry4BbgTvTwfRXSIKBtN59JIPofcAnI6IfoNA5s7oGMzMbmQ7HgFeldXV1RXd3d6WbYWZ2RJH0eER0jVTPU6SYmVlZHCRmZlYWB4mZmZXFQWJmZmWpisF2SVuAX43y8KOBl8ewOUeSar52qO7rr+Zrh+q+/txrf2NEdI50QFUESTkkdRfz1MJEVM3XDtV9/dV87VDd1z+aa3fXlpmZlcVBYmZmZXGQjGxxpRtQQdV87VDd11/N1w7Vff0lX7vHSMzMrCy+IzEzs7I4SMzMrCwOkkOQtEDSs5LWSbqu0u05nCRtkPSUpNWSJvyMl5KWSNosaZCK8KkAAAQkSURBVE1O2VGSfiTpufRnRyXbmJVhrv3PJW1Kf/+rJZ1byTZmRdIJkn4i6WlJayV9Oi2f8L/7Q1x7yb97j5EMQ1It8O/A2SRL+q4ELo2IqlgfXtIGoCsiquJLWZLOAl4D7oiIU9Oy/wW8EhE3pv+R6IiIz1eynVkY5tr/HHgtIv6ykm3LmqTjgOMiYpWkVuBx4HzgCib47/4Q134xJf7ufUcyvPnAuohYHxF7gaXAwgq3yTISEf+XZE2cXAuB29Pt20n+kk04w1x7VYiIlyJiVbq9A3gGmEYV/O4Pce0lc5AMbxqwMed9D6P8Qz5CBfBPkh6XdE2lG1Mhx0bES+n2r4FjK9mYCrhW0pNp19eE69rJJ2kGcBrwGFX2u8+7dijxd+8gseH8XkS8HXgv8Mm0+6NqpUtAV1M/8DeBNwFzgJeAv6psc7IlaTLwAPCZiHg1d99E/90XuPaSf/cOkuFtAk7IeT89LasKEbEp/bkZeJCkq6/a/CbtRx7sT95c4fYcNhHxm4joj4gB4BYm8O9fUj3JP6R3RcR30uKq+N0XuvbR/O4dJMNbCZwsaaakBpL15JdVuE2HhaRJ6eAbkiYB5wBrDn3UhLQMuDzdvhz4XgXbclgN/iOa+gMm6O9fkoBbgWci4n/n7Jrwv/vhrn00v3s/tXUI6WNvXwdqgSUR8ZUKN+mwkPRbJHchAHXA3RP92iXdA7yLZArt3wBfAr4L3AecSLIMwcURMeEGpYe59neRdG0EsAH4eM6YwYQh6feAfwGeAgbS4j8hGSuY0L/7Q1z7pZT4u3eQmJlZWdy1ZWZmZXGQmJlZWRwkZmZWFgeJmZmVxUFiZmZlcZCYjXOS3iXpB5Vuh9lwHCRmZlYWB4nZGJH0YUk/T9dwuFlSraTXJH0tXe/hnyV1pnXnSHo0nRjvwcGJ8SSdJOlhSb+QtErSm9LTT5Z0v6RfSror/Vay2bjgIDEbA5JmAZcAZ0TEHKAfuAyYBHRHxCnAT0m+NQ5wB/D5iPgdkm8WD5bfBdwUEW8Dfpdk0jxIZmb9DDAb+C3gjMwvyqxIdZVugNkE8R5gLrAyvVloJpnobwC4N63z98B3JLUB7RHx07T8duAf0vnNpkXEgwARsRsgPd/PI6Infb8amAH8a/aXZTYyB4nZ2BBwe0R84YBC6Yt59UY7J9GenO1+/HfXxhF3bZmNjX8GLpR0DAyt+f1Gkr9jF6Z1PgT8a0RsB7ZJOjMt/wjw03SVuh5J56fnaJTUclivwmwU/L8aszEQEU9L+jOSVSVrgH3AJ4HXgfnpvs0k4yiQTE3+rTQo1gMfS8s/Atws6Yb0HBcdxsswGxXP/muWIUmvRcTkSrfDLEvu2jIzs7L4jsTMzMriOxIzMyuLg8TMzMriIDEzs7I4SMzMrCwOEjMzK8v/B+OB7wAoe35rAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qDY-OczRNyFr"
      },
      "source": [
        "Visualization of the mean absolute error"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        },
        "id": "etCQ6U47NvnO",
        "outputId": "3b095070-6fff-4db1-aea7-8d0b787ba430"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "plt.plot(model.history.history['mean_absolute_error'])\n",
        "plt.title('mean absolute error')\n",
        "plt.ylabel('mean absolute error')\n",
        "plt.xlabel('epoch')\n",
        "plt.legend(['train', 'test'], loc='upper left')\n",
        "plt.show()"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5xdZX3v8c939twymZDbHhBCQiJEJSKGOlzkorSCBqnAqYqgWOTYppwjPXqsPWLrlR77Qk9r64VW0WLBA1KUUtMaS8EalMPFBIhIApgQgpmAkPttMvff+WOtyaxMJjNrhuzZO7O/79drv/Zaz3rW2r81O5nfrPWs53kUEZiZmY2kptwBmJnZ4cEJw8zMcnHCMDOzXJwwzMwsFycMMzPLxQnDzMxyccIwy0HSuZLaSnDckHTCoT6uWSk4YZgdJpxcrNycMMzsAJIKg9ZrR7n/qOrb4cEJw8adpPWS/lTS45L2SPoHSUdJ+pGkXZLulTQ9U/8MSQ9I2i7pF5LOzWy7StKT6X7rJP1RZtu5ktok/YmklyS9IOmqYeI66LEydf5M0ub0HN6XKX+7pNXpvhslfSyz7Q8lrZW0VdISSccc5POXSfqDzPoHJN2fLv80Lf6FpN2S3pOW/66klenP5gFJJw9zfq+RdE8ax9OSLs1s+0dJfy9pqaQ9wG+n5/hxSY8DeyTVSrpI0qr085ZJOjFzjAPqHywWO0xFhF9+jesLWA88BBwFzAJeAh4FTgEagf8EPpPWnQVsAd5O8gfO+el6S7r9QuB4QMCbgXbgt9Jt5wI9wHVAXXqMdmD6QeLKc6wvAQ3p9j3Aq9PtLwDnpMvTM/v9DrAZ+K10v68CP818ZgAnpMvLgD/IbPsAcP9QddP1U9Kf3elAAbgy/dk2DHFuk4ENwFVAbbrvZmBBuv0fgR3AWenPuTE91kpgNjAJeFV6zuenP8//BawF6jPf67765f535tehf/kKw8rlqxHxYkRsBH4GPBwRj0VEB3AXyS80gCuApRGxNCL6IuIeYAXJL38i4ocR8Uwk7gP+Azgn8zndwHUR0R0RS4HdwKuHCijHsQA+FRGd6fYfAv1/pXcDCyQdERHbIuLRtPx9wE0R8WhEdAKfAN4oae4YfmaDLQa+EREPR0RvRNwMdAJnDFH3d4H1EfHtiOiJiMeAO4F3Z+r8ICL+X/pz7kjLvhIRGyJiL/Ae4IcRcU9EdAN/RZJIzswcI1vfJhgnDCuXFzPLe4dYb06XjwPend4C2S5pO3A2cDSApAskPZTeZtlOkkiKmWNtiYiezHp75tj7yXGsbRGxJ7P+HNB/e+mdaf3nJN0n6Y1p+TFpPQAiYjfJFdKsoWIYpeOAPxn0s5mdiWlw3dMH1X0f8IpMnQ1D7JctG3wufen2WQepbxOM7zFapdsAfCci/nDwBkkNJH8l/z7JX8fdkv6F5JbSqOQ81nRJkzNJYw7wBEBELAcullQHXAPcQfLL+3mSX9b9nzMZmAlsHCKMPUBTZv0VQ9TJ2gB8PiI+n+MUNwD3RcT5w9QZaujqbNnzwOv6VySJ5Bw3HqS+TTC+wrBK93+Bd0h6m6SCpMa0MftYoJ6kXWAT0CPpAuCtY/ycvMf6nKR6SeeQ3Ob5Xrr+PklT01s1O4G+tP53gaskLUyT0l+S3H5bP8SxVwK/J6lJyeOzHxy0/UXglZn1bwJXSzpdicmSLpQ0ZYhj/xvwKknvl1SXvk7NNlrncAdwoaS3pInxT0hugT0wimPYYcwJwypaRGwALgb+jOSX+QbgT4GaiNgF/A+SX2TbgPcCS8b4OXmO9Zt02/PArcDVEfFUuu39wHpJO4GrSW73EBH3Ap8iuXp5gaRR/bKDhPE3QBdJYrg5/YyszwI3p7eULo2IFcAfAl9L41pL0lB+sPN7a/rZz6fn8gWSJJlLRDxN0qb0VZIG83cA74iIrrzHsMObInwFaWZmI/MVhpmZ5eKEYWZmuThhmJlZLk4YZmaWy4Tph1EsFmPu3LnlDsPM7LDyyCOPbI6Iljx1J0zCmDt3LitWrCh3GGZmhxVJz41cK+FbUmZmlosThpmZ5eKEYWZmuUyYNoyhdHd309bWRkdHx8iVD3ONjY0ce+yx1NXVlTsUM5ugJnTCaGtrY8qUKcydO5dkYM2JKSLYsmULbW1tzJs3r9zhmNkENaFvSXV0dDBz5swJnSwAJDFz5syquJIys/KZ0AkDmPDJol+1nKeZlc+ETxgj6e3r48WdHbR39Yxc2cysilV9wgjgxZ0d7OnsLcnxt2/fzt/93d+Ner+3v/3tbN++vQQRmZmNTdUnjIKEJHr6+kauPAYHSxg9PcNf0SxdupRp06aVJCYzs7GY0E9J5SGJ2hrR01uaiaSuvfZannnmGRYuXEhdXR2NjY1Mnz6dp556il/96ldccsklbNiwgY6ODj784Q+zePFiYGCok927d3PBBRdw9tln88ADDzBr1ix+8IMfMGnSpJLEa2Z2MFWTMD73r6tY/fzOIbft7e5FQGNdYVTHXHDMEXzmHa8dts7111/PE088wcqVK1m2bBkXXnghTzzxxL7HX2+66SZmzJjB3r17OfXUU3nnO9/JzJkz9zvGmjVr+O53v8s3v/lNLr30Uu68806uuOKKUcVqZvZyVU3CGI5I2jLGw2mnnbZfX4mvfOUr3HXXXQBs2LCBNWvWHJAw5s2bx8KFCwF4wxvewPr168cpWjOzASVNGJIWAV8GCsC3IuL6QduvBj4E9AK7gcURsVrSXOBJ4Om06kMRcfXLiWW4K4G2re3s6uzhxKOPeDkfkcvkyZP3LS9btox7772XBx98kKamJs4999wh+1I0NDTsWy4UCuzdu7fkcZqZDVayhCGpANwAnA+0AcslLYmI1Zlqt0XE19P6FwFfAhal256JiIWlii+rtpC0YUTEIe/PMGXKFHbt2jXkth07djB9+nSampp46qmneOihhw7pZ5uZHUqlvMI4DVgbEesAJN0OXAzsSxgRkW1UmMz43RnaT6GmhiDo7QtqC4c2YcycOZOzzjqLk046iUmTJnHUUUft27Zo0SK+/vWvc+KJJ/LqV7+aM84445B+tpnZoVTKhDEL2JBZbwNOH1xJ0oeAjwL1wO9kNs2T9BiwE/hkRPxsiH0XA4sB5syZM+ZA69Ik0dMX1I6u3TuX2267bcjyhoYGfvSjHw25rb+dolgs8sQTT+wr/9jHPnbI4zMzy6Ps/TAi4oaIOB74OPDJtPgFYE5EnEKSTG6TdEADQ0TcGBGtEdHa0pJrhsEh1dakCaO3NH0xzMwmglImjI3A7Mz6sWnZwdwOXAIQEZ0RsSVdfgR4BnhVieKktpD8GHr6ynJHzMzssFDKhLEcmC9pnqR64DJgSbaCpPmZ1QuBNWl5S9pojqRXAvOBdWMJImLkJDBwhXH4Jow852lm9nKUrA0jInokXQPcTfJY7U0RsUrSdcCKiFgCXCPpPKAb2AZcme7+JuA6Sd1AH3B1RGwdbQyNjY1s2bJlxCHOCzVClG54kFLrnw+jsbGx3KGY2QSmifKXaWtra6xYsWK/stHMuPfCjg4a62qY3lRfqhBLyjPumdlYSHokIlrz1J3QPb3r6upyz0D3J1/+Ga+Y2shNH3h9iaMyMzs8lf0pqUpRnNLA5t2d5Q7DzKxiOWGkis31bN7lhGFmdjBOGKmW5gY27+7y00ZmZgfhhJEqNjfQ1dvHzg5P1WpmNhQnjFRxSvJ0lNsxzMyG5oSRKjYnQ4i7HcPMbGhOGKl9CWN3V5kjMTOrTE4YqYGE4SsMM7OhOGGkZkyup0ZOGGZmB+OEkSrUiBmT650wzMwOwgkjo9jcwKZdbsMwMxuKE0ZGsdnDg5iZHYwTRkax2bekzMwOxgkjo/8Kw8ODmJkdyAkjozilgY7uPvZ09ZY7FDOziuOEkeHe3mZmB+eEkVFs9nhSZmYHU9KEIWmRpKclrZV07RDbr5b0S0krJd0vaUFm2yfS/Z6W9LZSxtnPvb3NzA6uZAlDUgG4AbgAWABcnk0Iqdsi4nURsRD4IvCldN8FwGXAa4FFwN+lxyuplilJwtjk8aTMzA5QyiuM04C1EbEuIrqA24GLsxUiYmdmdTLQ/3jSxcDtEdEZEc8Ca9PjldSMyektKbdhmJkdoLaEx54FbMistwGnD64k6UPAR4F64Hcy+z40aN9ZpQlzQF2hhulNdb4lZWY2hLI3ekfEDRFxPPBx4JOj2VfSYkkrJK3YtGnTIYnHvb3NzIZWyoSxEZidWT82LTuY24FLRrNvRNwYEa0R0drS0vIyw00U07m9zcxsf6VMGMuB+ZLmSaonacRekq0gaX5m9UJgTbq8BLhMUoOkecB84OcljHWf4hRfYZiZDaVkbRgR0SPpGuBuoADcFBGrJF0HrIiIJcA1ks4DuoFtwJXpvqsk3QGsBnqAD0XEuHS/LjbXu9HbzGwIpWz0JiKWAksHlX06s/zhYfb9PPD50kU3tGJzA3u6etnb1cuk+pI/yWtmdtgoe6N3pWlx5z0zsyE5YQxSnJL0xdjkhGFmth8njEE8AKGZ2dCcMAYZGE/Kj9aamWU5YQwy0yPWmpkNyQljkIbaAkc01jphmJkN4oQxBHfeMzM7kBPGEIrNDWze5TYMM7MsJ4whtHgAQjOzAwybMCQVJP3P8QqmUhSb690Pw8xskGETRjp+0+XjFEvFKDY3sKujh47ucRm+yszssJBnLKn/J+lrwD8Be/oLI+LRkkVVZsV0qtYte7qYNW1SmaMxM6sMeRLGwvT9ukxZMDA73oST7e3thGFmlhgxYUTEb49HIJWk6M57ZmYHGPEpKUlTJX2pfypUSX8taep4BFcuRY9Ya2Z2gDyP1d4E7AIuTV87gW+XMqhya5ni8aTMzAbL04ZxfES8M7P+OUkrSxVQJWisK9DcUMsmj1hrZrZPniuMvZLO7l+RdBawt3QhVYZic71vSZmZZeS5wrgauCXTbrFv7u2JrOje3mZm+xmxpzfw/oh4PXAycHJEnBIRj+c5uKRFkp6WtFbStUNs/6ik1ZIel/RjScdltvVKWpm+lozyvF62JGG4DcPMrF+ent5np8s7I2Jn3gOnyeYG4AJgAXC5pAWDqj0GtEbEycD3gS9mtu2NiIXp66K8n3uoFKf4lpSZWVaeW1KPpX/hf4/9e3r/8wj7nQasjYh1AJJuBy4GVmeO8ZNM/YeAK3LGXXLF5ga2t3fT3dtHXcFjNJqZ5flN2AhsIenZ/Y709bs59psFbMist6VlB/NB4EfZz037fTwk6ZKhdpC0uL9/yKZNm3KElF9/X4wtvi1lZgaMcIWR3lbaEhEfK2UQkq4AWoE3Z4qPi4iNkl4J/KekX0bEM9n9IuJG4EaA1tbWOJQxZTvvvWJq46E8tJnZYSlPG8ZZYzz2RmB2Zv3YtGw/ks4D/hy4KCL2NRpExMb0fR2wDDhljHGMScuUZHgQD3NuZpbI04axcoxtGMuB+ZLmkSSKy4D3ZitIOgX4BrAoIl7KlE8H2iOiU1KRJGllG8RLLjsAoZmZ5UsY2TaMfgEMmzAiokfSNcDdQAG4KSJWSboOWBERS4D/AzQD35ME8Ov0iagTgW9I6iO5Cro+IlYP+UElMnBLym0YZmaQb7Taq8Z68IhYCiwdVPbpzPJ5B9nvAeB1Y/3cQ2FyQy2T6gp+tNbMLJVntNpXpZ3qnkjXT5b0ydKHVn7ui2FmNiDPY7XfBD4BdAOkvbwvK2VQlcLDg5iZDciTMJoi4ueDynpKEUylKTY3sHmX2zDMzCBfwtgs6XiShm4kvQt4oaRRVQhfYZiZDcjzlNSHSDrHvUbSRuBZ4H0ljapCtDTXs7W9i57ePmo9PIiZVbk8T0mtA86TNBmoiYhdpQ+rMhSnNBABW9u7OHKKe3ubWXXL/WdzROyppmQB2c57bscwM/N9lmFkx5MyM6t2ThjDKDYn40k5YZiZ5eu41yTpU5K+ma7Pl5RnePPDXnGKrzDMzPrlucL4NtAJvDFd3wj875JFVEGmNNRSX1vj8aTMzMiXMI6PiC8y0NO7HVBJo6oQkmhpbvCItWZm5EsYXZImMdBx73iSK46qUGyu95wYZmbk67j3WeDfgdmSbiWZm2LMI9geborNDTy/o6PcYZiZlV2ejnv/IekR4AySW1EfjojNJY+sQhSbG3h8445yh2FmVnZ5npL6cURsiYgfRsS/RcRmST8ej+AqQXFKPVv3dNHXd0inDDczO+wc9ApDUiPQBBTTKVP7G7qPAGaNQ2wVodjcQG9fsK29i5lpRz4zs2o03C2pPwI+AhwDPJop3wl8rZRBVZLsVK1OGGZWzQ56SyoivhwR84CPRcS8zOv1EZErYUhaJOlpSWslXTvE9o9KWi3p8XRWv+My266UtCZ9XTmmszsEPDyImVkiz1NSOyT9/uDCiLhluJ0kFYAbgPOBNmC5pCURsTpT7TGgNSLaJf034IvAeyTNAD4DtJI8zvtIuu+2XGd1CLVM8fAgZmaQrx/GqZnXOSSP2V6UY7/TgLURsS4iuoDbgYuzFSLiJ2lHQICHgGPT5bcB90TE1jRJ3AMsyvGZh1z/FcYmd94zsyqX57HaP86uS5pG8st/JLOADZn1NuD0Yep/EPjRMPse0NAuaTGwGGDOnDk5Qhq9qZPqqCvIw4OYWdUby2i1e4B5hzIISVeQ3H76P6PZLyJujIjWiGhtaWk5lCFlY2PmZE/VamY24hWGpH8lHRaEJMEsAO7IceyNwOzM+rFp2eDjnwf8OfDmiOjM7HvuoH2X5fjMkihOqXfCMLOql6fR+68yyz3AcxHRlmO/5cB8SfNIEsBlwHuzFSSdAnwDWBQRL2U23Q38Zdr/A+CtwCdyfGZJFJt9hWFmlqcN476xHDgieiRdQ/LLvwDcFBGrJF0HrIiIJSS3oJqB70kC+HVEXBQRWyX9BUnSAbguIraOJY5DodjcwFMvVNXstGZmBxiup/cuBm5F7bcJiIg4YqSDR8RSYOmgsk9nls8bZt+bgJtG+ozxUGxuYMueTiKCNLGZmVWdgyaMiJgynoFUsmJzPd29wY693Uxrqi93OGZmZZGnDQNJryfpgwHw04h4vHQhVZ6WzFStThhmVq3yjFb7YeBW4Mj0daukPx5+r4lloPOe+2KYWfXKc4XxQeD0iNgDIOkLwIPAV0sZWCXxeFJmZvk67gnozaz3UiVzevcrNns8KTOzPFcY3wYelnQXSaK4GPiHkkZVYaY31VOokROGmVW1PP0wviRpGXB2WnRVRDxW0qgqTE2NmDG5ns1uwzCzKpZnaJDjgVUR8aik3wbOkfRsRGwvfXiVw729zaza5WnDuBPolXQC8HWS8aFuK2lUFajY7PGkzKy65UkYfRHRA/we8LWI+FPg6NKGVXlamhs8xLmZVbU8CaNb0uXA7wP/lpbVlS6kylSc0sCm3cnwIGZm1ShPwrgKeCPw+Yh4Nh199julDavyFJvr6erpY1dnT7lDMTMrixETRjoH98eAVZJeB2yMiC+UPLIKs6/znqdqNbMqlWdokAuBZ4CvAF8D1kq6oNSBVZqB3t5uxzCz6pSn495fA78dEWth32O2P2Rg/u2q4OFBzKza5WnD2NWfLFLrgKqbTag4xcODmFl1G24Cpd9LF1dIWkoyj3cA72ZgJryqMaOpHsltGGZWvYa7JfWOzPKLwJvT5U3ApJJFVKFqCzXMaKpnk9swzKxKDTfj3lUv9+CSFgFfJpnT+1sRcf2g7W8C/hY4GbgsIr6f2dYL/DJd/XVEXPRy43m5PDyImVWzPGNJNZLMifFaoLG/PCL+6wj7FYAbgPOBNmC5pCXpY7r9fg18gOSx3cH2RsTCkeIbT8UpHh7EzKpXnkbv7wCvAN4G3AccS75G79OAtRGxLiK6gNtJhkbfJyLWp9O99o0q6jLxFYaZVbM8CeOEiPgUsCcibgYuBE7Psd8sYENmvS0ty6tR0gpJD0m6ZKgKkhandVZs2rRpFIcem5mTGzzEuZlVrVxjSaXv2yWdBEwlmdu71I6LiFbgvcDfpv0/9hMRN0ZEa0S0trS0lDyg4pR69nb3ssfDg5hZFcqTMG6UNB34JLAEWA3kGRpkI8lQ6P2OTctyiYiN6fs6YBlwSt59S8Wd98ysmuUZS+pbEbEtIn4aEa+MiCMj4hs5jr0cmC9pnqR64DKShDMiSdMlNaTLReAskkRVVi1OGGZWxfJcYYxJOofGNcDdwJPAHRGxStJ1ki4CkHSqpDaSzoDfkLQq3f1Ekg6DvwB+Alw/6Omqsui/wtjkdgwzq0J5xpIas4hYCiwdVPbpzPJykltVg/d7AHhdKWMbCw8PYmbVrGRXGBPRzMm+JWVm1SvXFYakM4G52foRcUuJYqpY9bU1TJ1U54RhZlUpT0/v7wDHAyuB3rQ4gKpLGJDMvOe+GGZWjfJcYbQCC8KTWQPu7W1m1StPG8YTJEODGFCc4oRhZtUpzxVGEVgt6efAvt+UlTB6bDm0NDd4mlYzq0p5EsZnSx3E4aTYXM/uzh46untprCuUOxwzs3EzYsKIiPvGI5DDxUDnvU5mz2gqczRmZuNnxDYMSWdIWi5pt6QuSb2Sdo5HcJXI40mZWbXK0+j9NeByYA3J1Kx/QDIxUlUqTulPGG7HMLPqkqund0SsBQoR0RsR3wYWlTasylVs9vAgZlad8jR6t6ejza6U9EXgBap4SJF9t6R2OWGYWXXJ84v//Wm9a4A9JHNcvLOUQVWyxroCUxpqfYVhZlUnz1NSz0maBBwdEZ8bh5gqXtJ5z20YZlZd8jwl9Q6ScaT+PV1fKCnXREgTVbG5nk2+wjCzKpPnltRngdOA7QARsRKYV8KYKp7HkzKzapQnYXRHxI5BZVU9EGGxucGN3mZWdfIkjFWS3gsUJM2X9FXggRLHVdGKzQ3s7Oihs6d35MpmZhNEnoTxx8BrSQYe/C6wE/hInoNLWiTpaUlrJV07xPY3SXpUUo+kdw3adqWkNenryjyfN17mFpMhQf5+2TNljsTMbPzkeUqqHfjz9JWbpAJJj/DzgTZguaQlEbE6U+3XwAeAjw3adwbwGZK5OAJ4JN1322hiKJXfPfkYfrZmM3977xoAPnLeq8ockZlZ6eWZca8V+DMOnKL15BF2PQ1YGxHr0uPcDlwM7EsYEbE+3dY3aN+3AfdExNZ0+z0kvcu/O1K846FQI77wzuT0nTTMrFrk6el9K/CnwC+Bwb/YhzML2JBZbwNOfxn7zhpcSdJiYDHAnDlzRhHay+ekYWbVJk/C2BQRFdnvIiJuBG4EaG1tHfcnt5w0zKya5EkYn5H0LeDH7D/j3j+PsN9GkmFE+h2bluWxETh30L7Lcu47rpw0zKxa5EkYVwGvAeoYuCUVwEgJYzkwX9I8kgRwGfDenHHdDfylpOnp+luBT+Tcd9w5aZhZNciTME6NiFeP9sAR0SPpGpJf/gXgpohYJek6YEVELJF0KnAXMB14h6TPRcRrI2KrpL8gSToA1/U3gFcqJw0zm+jyJIwHJC0Y9DhsLhGxFFg6qOzTmeXlJLebhtr3JuCm0X5mOTlpmNlElidhnEEyF8azJG0YAiLHY7VVyUnDzCaqPAmjamfXGysnDTObiHLNhzEegUw0ThpmNtHkucKwMRqcNBrrClz95uPLHJWZ2dhU7dzc46U/aZy/4Ci+fO8adrR3lzskM7MxccIYB4Ua8eG3zGdvdy/fe2TDyDuYmVUgJ4xxctKsqbQeN51bHnyO3r6qnn/KzA5TThjj6Moz5/Lrre0se/qlcodiZjZqThjjaNFJr+CoIxr4xwfWlzsUM7NRc8IYR3WFGt53+nH8bM1mntm0u9zhmJmNihPGOLv8tDnUF2q4xVcZZnaYccIYZy1TGrjw5KP5/iNt7OrwI7ZmdvhwwiiDK8+cy56uXu58pK3coZiZ5eaEUQYLZ0/j9bOnccuDz9HnR2zN7DDhhFEmHzjzONZt3sPP1m4udyhmZrk4YZTJ2193NMXmem5247eZHSacMMqkobbAe0+bw0+efonntuwpdzhmZiNywiij951xHAWJWx70CPJmVvlKmjAkLZL0tKS1kq4dYnuDpH9Ktz8saW5aPlfSXkkr09fXSxlnuRx1RCOLTnoFd6zYwJ7OnnKHY2Y2rJIlDEkF4AbgAmABcLmkBYOqfRDYFhEnAH8DfCGz7ZmIWJi+ri5VnOX2gTPnsqujh7se21juUMzMhlXKK4zTgLURsS4iuoDbgYsH1bkYuDld/j7wFkkqYUwV5w3HTee1xxzBLQ+uJ8KP2JpZ5SplwpgFZCd/aEvLhqwTET3ADmBmum2epMck3SfpnBLGWVaSuPLMufzqxd08+MyWcodjZnZQldro/QIwJyJOAT4K3CbpiMGVJC2WtELSik2bNo17kIfKRa8/hulNdR7F1swqWikTxkZgdmb92LRsyDqSaoGpwJaI6IyILQAR8QjwDPCqwR8QETdGRGtEtLa0tJTgFMZHY12By06bw71PvkjbtvZyh2NmNqRSJozlwHxJ8yTVA5cBSwbVWQJcmS6/C/jPiAhJLWmjOZJeCcwH1pUw1rK74ozjAPjOQ37E1swqU8kSRtomcQ1wN/AkcEdErJJ0naSL0mr/AMyUtJbk1lP/o7dvAh6XtJKkMfzqiNhaqlgrwaxpk3jrglfwT8s30NHdW+5wzMwOUFvKg0fEUmDpoLJPZ5Y7gHcPsd+dwJ2ljK0SXXnmXP591W9YsvJ5Lj119sg7mJmNo0pt9K5KZ7xyBq8+agr/+IAfsTWzyuOEUUH6H7Fd/cJOVjy3rdzhmJntxwmjwlxyyjEc0VjrR2zNrOI4YVSYpvpa3nPqbP79id/wmx0d5Q7HzGwfJ4wK9P4z5tIXwa0P+xFbM6scThgVaM7MJt7ymiP57s9/TWePH7E1s8rghFGhrjxzLpt3d/EHN6/gOw+u59nNe/zklJmVVUn7YdjYnX1Ckf9+7vH8YOXzfOoHq4Ckc98584ucPb/IWccXmT65vsxRmlk10UT5q7W1tTVWrFhR7jAOuYhg/ZZ27l+ziZ+t2cyDz2xhV2cPEpx0zFTOnl/knBOKvGHudBpqC+UO18wOM5IeiYjWXHWdMA4vPaJ51Q8AAAqKSURBVL19/KJtB/ev2cz9azfx2K+309MXNNbVcNq8mbxpfpG3nHgU84qTyx2qmR0GnDCqyK6Obh5et5X7127mp2s2sW7THgBOOLKZ8048ivMXHMUps6dRU1NV81KZWU5OGFVsw9Z27n3yRe5Z/SIPP7uV3r6g2FzPW16TJI+z5xdprPOtKzNLOGEYADvau1n2q5e4Z/WL3Pf0JnZ19tBYV8M581s4f8FR/M5rjqTY3FDuMM2sjEaTMPyU1AQ2tamOixfO4uKFs+jq6ePhZ7dwz+oXuXd1cgUiwRvmTOeNx89k9owmjp0+idnTmzh6aiO1BT9xbWb78xVGFYoIVj2/c9+tqydf2Elf5p9BoUYcM62RY6c1MXtGkkSOTd9nz2iipbnBbSJmE4RvSdmodPX08cKOvbRt28uGre1s2NbOhq17advWzoZte9m0q3O/+vW1NUybVMfkhlqa6gtMbqhl8r73WpoaCjQ31NJUX0tzQ4Gm+lomN9QyramOGZPrmdZUx7RJ9dTX+irGrNx8S8pGpb62huNmTua4mUM/itvR3Zskk23ttG1tp23bXnbs7WZPVy/tnT3s7uxh8+4untvazp7OHto7e9nT1bPfVctQmtMkMr2pfl8y6V/OvvcvT2uqo7mhFslXN2bl4IRhI2qsK3DCkc2ccGRz7n0igo7uPvZ09bCns4ddHT3s2NvNtvYutrV3s21PF9vau9jePlD23JZ2trV3sauj56DHrSuIqZPqmZ4mkGlNyXKSVOppqi/QUFtDfW0NDbXZ5YGy/vWGuhoaCsl6fW0NBd9mMxuWE4aVhCQm1ReYVF8Y9ZNYPb19bN/bzfZMctne3s32vcn69vYutu1JEs2Gre083paUd/X0vayYa5RcbdUXavZ7r+tfz5Q11BZoqKuhMX1vqK2hsS5JUP2Jat/6fvUGyvYt19bQkNatrZGvoKxilTRhSFoEfBkoAN+KiOsHbW8AbgHeAGwB3hMR69NtnwA+CPQC/yMi7i5lrFY5ags1FJsbRpVoIoK93b3s7eqls6ePrp4+Onv66Ozp3bfcla537reevHf3Ju9dg98Hbevs6WNXRw+be7qSY3X3Zd6TOi9Hjdh3FdSfN7LNjIPbHGO/fTV0MsoktcEJq7ZGIBBCSj6/f1kAEkqP3b+9UFNDXUHUFWqoLYi6mvS9kJTX1tRQV1tDXY2oTevUHHAcUVPDvnJJ1PSXp3XSj0eD9k3iSuKsyWzfV1cc8HkM+uxsfcuvZAlDUgG4ATgfaAOWS1oSEasz1T4IbIuIEyRdBnwBeI+kBcBlwGuBY4B7Jb0qIjzWtw1JEk31SUN7OfX2xX6JqaO7l47uwclqIMHsK+veP5l1du//T324X2z9m/r6IklqQxx7256u/Y7dH1tPXyRJJyAIIqAvkrIJ8jzMiPqT0AGJZ9+2/RNNslNmnyH2h2x55rPY/3vcf1u2XAOxZWLor9f/Wf1xLDhmKl+9/JRD8eMYVin/d50GrI2IdQCSbgcuBrIJ42Lgs+ny94GvKflJXQzcHhGdwLOS1qbHe7CE8Zq9bIWagVtxE0XEQCLpiyQpdvf10dMbdPcmV2A9vUFPXx9dPcl7d7qtpzepS38iyhwn0veBBBX0RdDbl34m7JfIgoFjJNuSOn19AwkuBsWb3YdM3exxYt9ytmwghoHjpT8PBo43UOfA/cnUgwMTcGSuD/erlzkemRj337b/z2jOjElj/XpHpZQJYxawIbPeBpx+sDoR0SNpBzAzLX9o0L6zBn+ApMXAYoA5c+YcssDNbMC+2zyZv4EnMXESouV3WD8IHxE3RkRrRLS2tLSUOxwzswmtlAljIzA7s35sWjZkHUm1wFSSxu88+5qZ2TgqZcJYDsyXNE9SPUkj9pJBdZYAV6bL7wL+M5KbdUuAyyQ1SJoHzAd+XsJYzcxsBCVrw0jbJK4B7iZ5rPamiFgl6TpgRUQsAf4B+E7aqL2VJKmQ1ruDpIG8B/iQn5AyMysvjyVlZlbFRjOW1GHd6G1mZuPHCcPMzHJxwjAzs1wmTBuGpE3Acy/jEEVg8yEK53Djc69e1Xz+1XzuMHD+x0VEro5sEyZhvFySVuRt+JlofO7Vee5Q3edfzecOYzt/35IyM7NcnDDMzCwXJ4wBN5Y7gDLyuVevaj7/aj53GMP5uw3DzMxy8RWGmZnl4oRhZma5VH3CkLRI0tOS1kq6ttzxjDdJ6yX9UtJKSRN6MC5JN0l6SdITmbIZku6RtCZ9n17OGEvpIOf/WUkb0+9/paS3lzPGUpE0W9JPJK2WtErSh9PyCf/9D3Puo/7uq7oNI513/Fdk5h0HLh807/iEJmk90BoRE74Dk6Q3AbuBWyLipLTsi8DWiLg+/YNhekR8vJxxlspBzv+zwO6I+KtyxlZqko4Gjo6IRyVNAR4BLgE+wAT//oc590sZ5Xdf7VcY++Ydj4guoH/ecZuAIuKnJMPoZ10M3Jwu30zyH2lCOsj5V4WIeCEiHk2XdwFPkkz7POG//2HOfdSqPWEMNe/4mH6Qh7EA/kPSI+kc6dXmqIh4IV3+DXBUOYMpk2skPZ7esppwt2QGkzQXOAV4mCr7/gedO4zyu6/2hGFwdkT8FnAB8KH0tkVVSmd7rLZ7tH8PHA8sBF4A/rq84ZSWpGbgTuAjEbEzu22if/9DnPuov/tqTxhVP3d4RGxM318C7iK5TVdNXkzv8fbf632pzPGMq4h4MSJ6I6IP+CYT+PuXVEfyC/PWiPjntLgqvv+hzn0s3321J4w8845PWJImp41gSJoMvBV4Yvi9JpzsvPJXAj8oYyzjrv+XZeq/MEG/f0kimRL6yYj4UmbThP/+D3buY/nuq/opKYD0UbK/ZWDe8c+XOaRxI+mVJFcVkMzvfttEPn9J3wXOJRnW+UXgM8C/AHcAc0iGx780IiZkw/BBzv9cklsSAawH/ihzT3/CkHQ28DPgl0BfWvxnJPfyJ/T3P8y5X84ov/uqTxhmZpZPtd+SMjOznJwwzMwsFycMMzPLxQnDzMxyccIwM7NcnDDMKoCkcyX9W7njMBuOE4aZmeXihGE2CpKukPTzdP6Ab0gqSNot6W/SuQZ+LKklrbtQ0kPp4G539Q/uJukESfdK+oWkRyUdnx6+WdL3JT0l6da0h65ZxXDCMMtJ0onAe4CzImIh0Au8D5gMrIiI1wL3kfSgBrgF+HhEnEzSy7a//Fbghoh4PXAmycBvkIwi+hFgAfBK4KySn5TZKNSWOwCzw8hbgDcAy9M//ieRDFbXB/xTWuf/Av8saSowLSLuS8tvBr6Xjt01KyLuAoiIDoD0eD+PiLZ0fSUwF7i/9Kdllo8Thll+Am6OiE/sVyh9alC9sY6305lZ7sX/P63C+JaUWX4/Bt4l6UjYNx/0cST/j96V1nkvcH9E7AC2STonLX8/cF8641mbpEvSYzRIahrXszAbI/8FY5ZTRKyW9EmSGQprgG7gQ8Ae4LR020sk7RyQDJf99TQhrAOuSsvfD3xD0nXpMd49jqdhNmYerdbsZZK0OyKayx2HWan5lpSZmeXiKwwzM8vFVxhmZpaLE4aZmeXihGFmZrk4YZiZWS5OGGZmlsv/B9hGvDeH1QHuAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bg41axVhUvo0"
      },
      "source": [
        "Processing\n",
        "--------"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FOifYwxfJ9JZ"
      },
      "source": [
        "VDSR Model\n",
        "\n",
        "The landscape of the test area (Austria) is formed by fluvial and glacial erosion. Thus the model may not fit to other landscapes.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YCUO3kpYJuDO"
      },
      "source": [
        "**Skip the next step when training was applied before!**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O6Jb-0eIJvKJ"
      },
      "source": [
        "from keras import Sequential\n",
        "from keras.layers import Conv2D\n",
        "from keras.models import Sequential, Model\n",
        "from keras.layers import Dense, Activation, Input, add\n",
        "from keras import optimizers\n",
        "from keras import losses\n",
        "from keras.optimizers import Adam\n",
        "import tensorflow as tf\n",
        "\n",
        "IMG_SIZE = (300, 300, 1)\n",
        "\n",
        "model = Sequential()\n",
        "\n",
        "input_img = Input(shape=IMG_SIZE)\n",
        "\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(input_img)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "model = Conv2D(32, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "model = Activation('relu')(model)\n",
        "\n",
        "model = Conv2D(1, (9, 9), padding='same', kernel_initializer='he_normal')(model)\n",
        "\n",
        "\n",
        "res_img = model\n",
        "\n",
        "output_img = add([res_img, input_img])\n",
        "\n",
        "model = Model(input_img, output_img)\n",
        "\n",
        "adam = Adam(lr=0.0001)\n",
        "model.compile(adam, loss='mse', metrics=[tf.keras.metrics.MeanAbsoluteError()])\n",
        "\n",
        "# adapt path\n",
        "model.load_weights('/content/drive/My Drive/vdsr4dem/weights/vdsr4dem_scale_3.h5')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "saMMvLSAKy3X"
      },
      "source": [
        "Read all GeoTiffs from folder and process "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3M4onRoQas4Y"
      },
      "source": [
        "Skip the next step if rasterio is allready installed"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JodA0kVbarew"
      },
      "source": [
        "!pip install rasterio"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J9M0T6yiMr92",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "b0c12c3f-e70e-47ef-d42b-4c555823a588"
      },
      "source": [
        "import os\n",
        "import sys\n",
        "import numpy as np\n",
        "import rasterio\n",
        "from tensorflow.python.data.experimental import AUTOTUNE\n",
        "from rasterio.plot import reshape_as_image, reshape_as_raster, show\n",
        "\n",
        "# adapt paths\n",
        "path = '/content/drive/My Drive/vdsr4dem/lr_test'\n",
        "path_out = '/content/drive/My Drive/vdsr4dem/sr_images'\n",
        "\n",
        "for file in os.listdir(path):\n",
        "    current = os.path.join(path, file)\n",
        "    print(current)\n",
        "    x, y = [], []\n",
        "    with rasterio.open(current) as src:\n",
        "        array = src.read(\n",
        "            out_shape=(\n",
        "                src.count,\n",
        "                int(src.height),\n",
        "                int(src.width)\n",
        "            )\n",
        "        )\n",
        "        profile = src.profile                \n",
        "    out_profile = src.profile.copy()\n",
        "    max_val = np.amax(array)\n",
        "    min_val = np.amin(array)\n",
        "    array = reshape_as_image(array)\n",
        "    lr = (array - min_val) / (max_val - min_val)\n",
        "    x = np.expand_dims(lr, axis=0)\n",
        "    out_array = model.predict(x)\n",
        "    srA = out_array\n",
        "    srA = srA * (max_val - min_val) + min_val\n",
        "    srA = np.squeeze(srA, axis=0)\n",
        "    #print(srA.shape)\n",
        "    print (\"sr :\", srA.dtype, \" shape \", srA.shape)\n",
        "    current_out = os.path.join(path_out, file)\n",
        "    with rasterio.open(current_out, 'w', **out_profile) as dst:\n",
        "        dst.write(reshape_as_raster(srA))\n",
        "print(\"finished\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_iHcBe9cMdfY"
      },
      "source": [
        "Statistics with test dataset\n",
        "--------------------------------------\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Mafz4IzJa1IJ"
      },
      "source": [
        "Skip the next step if rasterio is allready installed"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QHgiUAh7a6IK",
        "outputId": "27355845-e36b-4261-fb2f-b51497073f16"
      },
      "source": [
        "!pip install rasterio"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
           
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FzFNr5E3AMIl",
        "outputId": "d58c6f4a-a5b2-4c2b-b26f-fff8053dcaa6"
      },
      "source": [
        "import os\n",
        "import sys\n",
        "import numpy as np\n",
        "from sklearn.metrics import mean_absolute_error\n",
        "from sklearn.metrics import mean_squared_error\n",
        "from skimage.metrics import peak_signal_noise_ratio\n",
        "from skimage.metrics import structural_similarity\n",
        "from keras.preprocessing.image import save_img\n",
        "import rasterio\n",
        "from tensorflow.python.data.experimental import AUTOTUNE\n",
        "from rasterio.plot import reshape_as_image, reshape_as_raster, show\n",
        "from rasterio import Affine\n",
        "from rasterio.enums import Resampling\n",
        "from numpy import savetxt\n",
        "\n",
        "# adapt paths\n",
        "path_lr = '/content/drive/My Drive/vdsr4dem/lr_test'\n",
        "path_hr = '/content/drive/My Drive/vdsr4dem/hr_test'\n",
        "\n",
        "meanAbsoluteError = []\n",
        "meanAbsoluteErrorCubic = []\n",
        "\n",
        "countMAE = 0\n",
        "countRMSE = 0\n",
        "countPSNR = 0\n",
        "countSSI = 0\n",
        "countMAESlice = 0\n",
        "countRMSESlice = 0\n",
        "countPSNRSlice = 0\n",
        "countSSISlice = 0\n",
        "totalMAE_lr = 0\n",
        "totalMAE_sr = 0\n",
        "totalRMSE_lr = 0\n",
        "totalRMSE_sr = 0\n",
        "totalMAX_lr = 0\n",
        "totalMAX_sr = 0\n",
        "fileCounter = 0 \n",
        "\n",
        "# read low resolution images\n",
        "for file in os.listdir(path_lr):\n",
        "    current = os.path.join(path_lr, file)\n",
        "    fileCounter += 1;\n",
        "    x, y = [], []\n",
        "    with rasterio.open(current) as src:\n",
        "        arrayLR = src.read(\n",
        "            out_shape=(\n",
        "                src.count,\n",
        "                int(src.height),\n",
        "                int(src.width)\n",
        "            )\n",
        "        )        \n",
        "    max_val_lr = np.amax(arrayLR)\n",
        "    min_val_lr = np.amin(arrayLR)\n",
        "    arrayLR = reshape_as_image(arrayLR)\n",
        "    nsamples, nx, ny = arrayLR.shape\n",
        "    arrayLR = arrayLR.reshape((nsamples,nx*ny))\n",
        "        \n",
        "    # process super resolution image\n",
        "    # normalize\n",
        "    lr4sr = (arrayLR - min_val_lr) / (max_val_lr - min_val_lr)\n",
        "    lr4sr = np.expand_dims(lr4sr, axis=0)\n",
        "    arraySR = model.predict(lr4sr)\n",
        "    # recalculate real world values for sr image\n",
        "    arraySR = arraySR * (max_val_lr - min_val_lr) + min_val_lr\n",
        "    arraySR = np.squeeze(arraySR, axis=0)\n",
        "    nsamples, nx, ny = arraySR.shape\n",
        "    arraySR = arraySR.reshape((nsamples,nx*ny))\n",
        "\n",
        "    # maximal value of sr image    \n",
        "    max_val_sr = np.amax(arraySR)\n",
        "\n",
        "    # load hr image to compare images    \n",
        "    currentHr = os.path.join(path_hr, file)\n",
        "    with rasterio.open(currentHr) as src:\n",
        "        arrayHR = src.read(\n",
        "            out_shape=(\n",
        "                src.count,\n",
        "                int(src.height),\n",
        "                int(src.width)\n",
        "            )\n",
        "        )\n",
        " \n",
        "    arrayHR = reshape_as_image(arrayHR)\n",
        "\n",
        "    nsamples, nx, ny = arrayHR.shape\n",
        "    arrayHR = arrayHR.reshape((nsamples,nx*ny))\n",
        "\n",
        "    # slice 15 pixels from border because border effects\n",
        "    slice_val = 15\n",
        "    HrSlice = arrayHR[slice_val:-slice_val, slice_val:-slice_val]\n",
        "    SrSlice = arraySR[slice_val:-slice_val, slice_val:-slice_val]\n",
        "    LrSlice = arrayLR[slice_val:-slice_val, slice_val:-slice_val]\n",
        "    \n",
        "    # calculate the metrics for sr image\n",
        "    MAE_sr = mean_absolute_error(arrayHR, arraySR)\n",
        "    RMSE_sr = mean_squared_error(arrayHR, arraySR, squared=False)\n",
        "    SSI_sr = structural_similarity(arrayHR, arraySR)\n",
        "    MAE_sr_slice = mean_absolute_error(HrSlice, SrSlice)\n",
        "    RMSE_sr_slice = mean_squared_error(HrSlice, SrSlice, squared=False)\n",
        "    SSI_sr_slice = structural_similarity(HrSlice, SrSlice)\n",
        "    MAX_sr = np.amax(abs(arrayHR - arraySR))\n",
        "    MAX_sr_slice = np.amax(abs(HrSlice - SrSlice))\n",
        "    # adapt data_range to the maximal overall value  of the dataset\n",
        "    PSNR_sr = peak_signal_noise_ratio(arrayHR, arraySR, data_range=3150)\n",
        "    PSNR_sr_slice = peak_signal_noise_ratio(HrSlice, SrSlice, data_range=3150)\n",
        "    \n",
        "    # calculate the metrics for lr image\n",
        "    MAE_lr = mean_absolute_error(arrayHR, arrayLR)\n",
        "    RMSE_lr = mean_squared_error(arrayHR, arrayLR, squared=False)\n",
        "    SSI_lr = structural_similarity(arrayHR, arrayLR)\n",
        "    MAE_lr_slice = mean_absolute_error(HrSlice, LrSlice)\n",
        "    RMSE_lr_slice = mean_squared_error(HrSlice, LrSlice, squared=False)\n",
        "    SSI_lr_slice = structural_similarity(HrSlice, LrSlice)\n",
        "    MAX_lr = np.amax(abs(arrayHR - arrayLR))\n",
        "    MAX_lr_slice = np.amax(abs(HrSlice - LrSlice))\n",
        "    # adapt data_range to the maximal overall value  of the dataset\n",
        "    PSNR_lr = peak_signal_noise_ratio(arrayHR, arrayLR, data_range=3150)\n",
        "    PSNR_lr_slice = peak_signal_noise_ratio(HrSlice, LrSlice, data_range=3150)\n",
        "    \n",
        "    # MAE count improved images\n",
        "    if MAE_sr < MAE_lr:\n",
        "        countMAE += 1\n",
        "    if MAE_sr_slice < MAE_lr_slice:\n",
        "        countMAESlice += 1    \n",
        "    # RMSE count improved images\n",
        "    if RMSE_sr < RMSE_lr:\n",
        "        countRMSE += 1\n",
        "    if RMSE_sr_slice < RMSE_lr_slice:\n",
        "        countRMSESlice += 1 \n",
        "    # PSNR count imporved images, higher values are better\n",
        "    if PSNR_sr > PSNR_lr:\n",
        "        countPSNR += 1\n",
        "    if PSNR_sr_slice > PSNR_lr_slice:\n",
        "        countPSNRSlice += 1    \n",
        "    # SSI count imporved images, higher values are better\n",
        "    if SSI_sr > SSI_lr:\n",
        "        countSSI += 1\n",
        "    if SSI_sr_slice > SSI_lr_slice:\n",
        "        countSSISlice += 1\n",
        "\n",
        "    totalMAE_lr = totalMAE_lr + MAE_lr_slice\n",
        "    totalMAE_sr = totalMAE_sr + MAE_sr_slice\n",
        "    totalRMSE_lr = totalRMSE_lr + RMSE_lr_slice\n",
        "    totalRMSE_sr = totalRMSE_sr + RMSE_sr_slice        \n",
        "    totalMAX_lr = totalMAX_lr + MAX_lr_slice\n",
        "    totalMAX_sr = totalMAX_sr + MAX_sr_slice \n",
        "\n",
        "    # print metrics for every tile if desired\n",
        "    #print(file ,' SR MaxError: ', MAX_sr_slice  ,' LR MaxError: ', MAX_lr_slice, ' Difference MaxError LR - SR: ', MAX_lr_slice - MAX_sr_slice)\n",
        "    #print(file ,' SR MAE: ', MAE_sr_slice ,' LR MAE: ', MAE_lr_slice, ' Difference LR - SR: ', MAE_lr_slice - MAE_sr_slice)\n",
        "    #print(file ,' SR RMSE: ', RMSE_sr_slice ,' LR RMSE: ', RMSE_lr_slice, ' Difference LR - SR: ', RMSE_lr_slice - RMSE_sr_slice)\n",
        "    #print(file ,' SR PSNR: ', PSNR_sr_slice ,' LR PSNR: ', PSNR_lr_slice, ' Difference LR - SR: ', PSNR_lr_slice - PSNR_sr_slice)\n",
        "\n",
        "print('Better results of',fileCounter,'test areas for MAE        ',countMAE)\n",
        "print('Better results of',fileCounter,'test areas for MAE sliced ',countMAESlice)\n",
        "print('Better results of',fileCounter,'test areas for RMSE       ',countRMSE)\n",
        "print('Better results of',fileCounter,'test areas for RMSE sliced',countRMSESlice)\n",
        "print('Better results of',fileCounter,'test areas for PNSR       ',countPSNR)\n",
        "print('Better results of',fileCounter,'test areas for PNSR sliced',countPSNRSlice)\n",
        "print('Better results of',fileCounter,'test areas for SSI        ',countSSI)\n",
        "print('Better results of',fileCounter,'test areas for SSI sliced ',countSSISlice)\n",
        "\n",
        "print('Mean MAE LR', totalMAE_lr/fileCounter);\n",
        "print('Mean MAE SR', totalMAE_sr/fileCounter);\n",
        "print('Difference LR - SR MAE ', totalMAE_lr/fileCounter - totalMAE_sr/fileCounter)\n",
        "print('Mean RMSE LR', totalRMSE_lr/fileCounter);\n",
        "print('Mean RMSE SR', totalRMSE_sr/fileCounter);\n",
        "print('Difference SR - LR RMSE ', totalRMSE_lr/fileCounter - totalRMSE_sr/fileCounter)\n",
        "print('MAX LR', totalMAX_lr/fileCounter);\n",
        "print('MAX SR', totalMAX_sr/fileCounter);\n",
        "print('Difference LR - SR MAX ', totalMAX_lr/fileCounter - totalMAX_sr/fileCounter)"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Better results of 30 test areas for MAE         30\n",
            "Better results of 30 test areas for MAE sliced  30\n",
            "Better results of 30 test areas for RMSE        30\n",
            "Better results of 30 test areas for RMSE sliced 30\n",
            "Better results of 30 test areas for PNSR        30\n",
            "Better results of 30 test areas for PNSR sliced 30\n",
            "Better results of 30 test areas for SSI         27\n",
            "Better results of 30 test areas for SSI sliced  30\n",
            "Mean MAE LR 0.855532372991244\n",
            "Mean MAE SR 0.8161334057648977\n",
            "Difference LR - SR MAE  0.03939896722634628\n",
            "Mean RMSE LR 1.480201146999995\n",
            "Mean RMSE SR 1.3770212570826212\n",
            "Difference SR - LR RMSE  0.1031798899173737\n",
            "MAX LR 24.78179931640625\n",
            "MAX SR 23.896238199869792\n",
            "Difference LR - SR MAX  0.8855611165364579\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}