notebooks/15_rnn_weather_forecasting.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "grBYZTvSVnpV"
      },
      "source": [
        "# Weather forecasting with RNNs\n",
        "\n",
        "This notebook is based on one of the Keras examples, where we will try to predict the weather using data recorded at the Weather Station of the Max Planck Institute for Biogeochemistry in Jena, Germany.\n",
        "\n",
        "Updated version of various additional data are available at www.bgc-jena.mpg.de/wetter -- the particular data we are looking at here, was recorded from 2009 to 2016."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "VH-vHuaxVtQ7"
      },
      "outputs": [],
      "source": [
        "import keras\n",
        "import tensorflow as tf\n",
        "import os\n",
        "import numpy as np\n",
        "from matplotlib import pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Uumj5I1zWATW"
      },
      "source": [
        "## Download the data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "c4OjAX7IWBYu"
      },
      "outputs": [],
      "source": [
        "!wget https://s3.amazonaws.com/keras-datasets/jena_climate_2009_2016.csv.zip\n",
        "!unzip jena_climate_2009_2016.csv.zip"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J2B8LZu2WHNo"
      },
      "source": [
        "Have a look at the contents. There are 14 features in total."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PhqxjSUUWIjE"
      },
      "outputs": [],
      "source": [
        "\n",
        "fname = os.path.join(\"jena_climate_2009_2016.csv\")\n",
        "with open(fname) as f:\n",
        "    data = f.read()\n",
        "    lines = data.split(\"\\n\")\n",
        "    header = lines[0].split(\",\")\n",
        "    lines = lines[1:]\n",
        "    print(header)\n",
        "    print(len(lines))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BlFoNpv4WbAt"
      },
      "source": [
        "Convert the data to numpy arrays:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "VwSoTmXoWcbb"
      },
      "outputs": [],
      "source": [
        "temperature = np.zeros((len(lines),))\n",
        "raw_data = np.zeros((len(lines), len(header) - 1))\n",
        "for i, line in enumerate(lines):\n",
        "    values = [float(x) for x in line.split(\",\")[1:]]\n",
        "    temperature[i] = values[1]\n",
        "    raw_data[i, :] = values[:]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kVem4AnZWiYA"
      },
      "source": [
        "What does the temperature look like?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Stxsq0fgWjsD"
      },
      "outputs": [],
      "source": [
        "plt.plot(range(len(temperature)), temperature)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "78X6Ye0RWtDz"
      },
      "source": [
        "A lot of data points here, so make a plot that focusses on the first 10 days. There is a measurement every 10 minutes, so we get 24 × 6 = 144 data points per day."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QjOM5bSAWvaF"
      },
      "outputs": [],
      "source": [
        "plt.plot(range(1440), temperature[:1440])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NaexxNJFW1dR"
      },
      "source": [
        "For our our experiments, we use the first 50% of the data for training, the following 25% for validation, and the last 25% for testing."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "v5-CHWvgW5FT"
      },
      "outputs": [],
      "source": [
        "num_train_samples = int(0.5 * len(raw_data))\n",
        "num_val_samples = int(0.25 * len(raw_data))\n",
        "num_test_samples = len(raw_data) - num_train_samples - num_val_samples\n",
        "print(\"num_train_samples:\", num_train_samples)\n",
        "print(\"num_val_samples:\", num_val_samples)\n",
        "print(\"num_test_samples:\", num_test_samples)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7P3amP1dXC1T"
      },
      "source": [
        "To prepare the data, we normalise it by subtracting the mean and dividing by the standard deviation."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XDbqingVXEIS"
      },
      "outputs": [],
      "source": [
        "mean = raw_data[:num_train_samples].mean(axis=0)\n",
        "raw_data -= mean\n",
        "std = raw_data[:num_train_samples].std(axis=0)\n",
        "raw_data /= std"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6pvmYXplXnxL"
      },
      "source": [
        "## Create the datasets to use\n",
        "\n",
        "We’ll use `timeseries_dataset_from_array()` to instantiate three datasets: one for training, one for validation, and one for testing.\n",
        "We’ll use the following parameter values:\n",
        "\n",
        " - `sampling_rate = 6`: Observations will be sampled at one data point per hour:\n",
        "we will only keep one data point out of 6.\n",
        " - `sequence_length = 120`: Observations will go back 5 days (120 hours).\n",
        " - `delay = sampling_rate * (sequence_length + 24 - 1)`: The target for a sequence will be the temperature 24 hours after the end of the sequence."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OfO5LNpFXq0j"
      },
      "outputs": [],
      "source": [
        "sampling_rate = 6\n",
        "sequence_length = 120\n",
        "delay = sampling_rate * (sequence_length + 24 - 1)\n",
        "batch_size = 256\n",
        "\n",
        "train_dataset = keras.utils.timeseries_dataset_from_array(\n",
        "    raw_data[:-delay],\n",
        "    targets=temperature[delay:],\n",
        "    sampling_rate=sampling_rate,\n",
        "    sequence_length=sequence_length,\n",
        "    shuffle=True,\n",
        "    batch_size=batch_size,\n",
        "    start_index=0,\n",
        "    end_index=num_train_samples - (num_train_samples % batch_size)\n",
        ")\n",
        "\n",
        "val_dataset = keras.utils.timeseries_dataset_from_array(\n",
        "    raw_data[:-delay],\n",
        "    targets=temperature[delay:],\n",
        "    sampling_rate=sampling_rate,\n",
        "    sequence_length=sequence_length,\n",
        "    shuffle=True,\n",
        "    batch_size=batch_size,\n",
        "    start_index=num_train_samples,\n",
        "    end_index=num_train_samples + num_val_samples\n",
        ")\n",
        "\n",
        "test_dataset = keras.utils.timeseries_dataset_from_array(\n",
        "    raw_data[:-delay],\n",
        "    targets=temperature[delay:],\n",
        "    sampling_rate=sampling_rate,\n",
        "    sequence_length=sequence_length,\n",
        "    shuffle=True,\n",
        "    batch_size=batch_size,\n",
        "    start_index=num_train_samples + num_val_samples)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "W09AWqMJYWpW"
      },
      "source": [
        "Try out the datasets:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9AYPdRpiYbXI"
      },
      "outputs": [],
      "source": [
        "for samples, targets in train_dataset:\n",
        "    print(\"samples shape:\", samples.shape)\n",
        "    print(\"targets shape:\", targets.shape)\n",
        "    break"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NiSKOvqUYjsn"
      },
      "source": [
        "## A simplistic model\n",
        "\n",
        "First, for the simplest approach -- 100% autocorrelation.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XzRRg-0gY0i6"
      },
      "outputs": [],
      "source": [
        "def evaluate_naive_method(dataset):\n",
        "    total_abs_err = 0.\n",
        "    samples_seen = 0\n",
        "    for samples, targets in dataset:\n",
        "        preds = samples[:, -1, 1] * std[1] + mean[1]\n",
        "        total_abs_err += np.sum(np.abs(preds - targets))\n",
        "        samples_seen += samples.shape[0]\n",
        "\n",
        "        return total_abs_err / samples_seen\n",
        "\n",
        "\n",
        "print(f\"Validation MAE: {evaluate_naive_method(val_dataset):.2f}\")\n",
        "print(f\"Test MAE: {evaluate_naive_method(test_dataset):.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WcXcyc0RY9_7"
      },
      "source": [
        "The most basic machine learning model -- a single-layer dense model."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NP3NH3wnZD8a"
      },
      "outputs": [],
      "source": [
        "inputs = keras.Input(shape=(sequence_length, raw_data.shape[-1]))\n",
        "print('inputs.shape:', inputs.shape)\n",
        "x = keras.layers.Reshape((sequence_length, raw_data.shape[-1]))(inputs)\n",
        "x = keras.layers.Dense(16, activation=\"relu\")(x)\n",
        "outputs = keras.layers.Dense(1)(x)\n",
        "\n",
        "model = keras.Model(inputs, outputs)\n",
        "\n",
        "callbacks = [\n",
        "    keras.callbacks.ModelCheckpoint(\"jena_dense.keras\",\n",
        "    save_best_only=True)\n",
        "]\n",
        "\n",
        "model.compile(optimizer=\"rmsprop\", loss=\"mse\", metrics=[\"mae\"])\n",
        "\n",
        "history = model.fit(\n",
        "    train_dataset,\n",
        "    epochs=10,\n",
        "    validation_data=val_dataset,\n",
        "    callbacks=callbacks\n",
        ")\n",
        "\n",
        "model = keras.models.load_model(\"jena_dense.keras\")\n",
        "\n",
        "print(f\"Test MAE: {model.evaluate(test_dataset)[1]:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fuT0mIQAckXC"
      },
      "source": [
        "Define a utulity function for plotting the learning curve"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jcdvI6gucl3-"
      },
      "outputs": [],
      "source": [
        "def plot_loss_curve(history):\n",
        "    loss = history.history[\"mae\"]\n",
        "    val_loss = history.history[\"val_mae\"]\n",
        "    epochs = range(1, len(loss) + 1)\n",
        "    plt.figure()\n",
        "    plt.plot(epochs, loss, \"bo\", label=\"Training MAE\")\n",
        "    plt.plot(epochs, val_loss, \"b\", label=\"Validation MAE\")\n",
        "    plt.title(\"Training and validation MAE\")\n",
        "    plt.legend()\n",
        "    plt.show()\n",
        "\n",
        "plot_loss_curve(history)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ztYxytgtc9ez"
      },
      "source": [
        "## A convolutional model\n",
        "\n",
        "Let's see how a straight-up convolutional model performs."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cyOtn5_rrSM4"
      },
      "source": [
        "### <span style=\"color: red; font-weight: bold;\">Exercise:<span>\n",
        "\n",
        "Build a model with two convolutional layers, each followed by a max-pooling layer. The convolutions can have 24 filters with kernel size 8."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iLIsurxLc_As"
      },
      "outputs": [],
      "source": [
        "inputs = keras.Input(shape=(sequence_length, raw_data.shape[-1]))\n",
        "x = ... # TODO\n",
        "x = keras.layers.GlobalAveragePooling1D()(x)\n",
        "outputs = keras.layers.Dense(1)(x)\n",
        "\n",
        "model = keras.Model(inputs, outputs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qdnsmv3Jrvb5"
      },
      "source": [
        "Now train it, and plot the learning curves:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GESCi0XNr1K0"
      },
      "outputs": [],
      "source": [
        "callbacks = [\n",
        "    keras.callbacks.ModelCheckpoint(\"jena_conv.keras\",\n",
        "    save_best_only=True)\n",
        "]\n",
        "\n",
        "model.compile(optimizer=\"rmsprop\", loss=\"mse\", metrics=[\"mae\"])\n",
        "history = model.fit(\n",
        "    train_dataset,\n",
        "    epochs=10,\n",
        "    validation_data=val_dataset,\n",
        "    callbacks=callbacks\n",
        ")\n",
        "model = keras.models.load_model(\"jena_conv.keras\")\n",
        "\n",
        "print(f\"Test MAE: {model.evaluate(test_dataset)[1]:.2f}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pohHaH8levai"
      },
      "outputs": [],
      "source": [
        "plot_loss_curve(history)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cow8IcHwiJO4"
      },
      "source": [
        "### <span style=\"color: red; font-weight: bold;\">Exercise:<span>\n",
        "\n",
        "Implement a simple [WaveNet](https://deepmind.google/discover/blog/wavenet-a-generative-model-for-raw-audio/) model -- that is, a fully convolutional network with strides larger than one, and padding set to `\"causal\"`, so that the layers can only look backwards in time.\n",
        "\n",
        "You can use the textbook to get a simple solution -- but try also to make a more complicated version.\n",
        "\n",
        "Then train and evaluate it."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "K-J3YY-1tCp6"
      },
      "outputs": [],
      "source": [
        "# model goes here\n",
        "\n",
        "wavenet = keras.Model(inputs, outputs)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QKWmcU9EtLIY"
      },
      "outputs": [],
      "source": [
        "wavenet.compile(optimizer=\"rmsprop\", loss=\"mse\", metrics=[\"mae\"])\n",
        "history = wavenet.fit(\n",
        "    train_dataset,\n",
        "    epochs=10,\n",
        "    validation_data=val_dataset,\n",
        "    callbacks=callbacks\n",
        ")\n",
        "\n",
        "print(f\"Test MAE: {wavenet.evaluate(test_dataset)[1]:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gl-T6HWNe597"
      },
      "source": [
        "# Compare to a recurrent model: LSTM"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZH9QRJLDe8KP"
      },
      "outputs": [],
      "source": [
        "inputs = keras.Input(shape=(sequence_length, raw_data.shape[-1]))\n",
        "x = keras.layers.LSTM(16)(inputs)\n",
        "outputs = keras.layers.Dense(1)(x)\n",
        "model = keras.Model(inputs, outputs)\n",
        "callbacks = [\n",
        "    keras.callbacks.ModelCheckpoint(\"jena_lstm.keras\",\n",
        "    save_best_only=True)\n",
        "]\n",
        "\n",
        "model.compile(optimizer=\"rmsprop\", loss=\"mse\", metrics=[\"mae\"])\n",
        "history = model.fit(\n",
        "    train_dataset,\n",
        "    epochs=10,\n",
        "    validation_data=val_dataset,\n",
        "    callbacks=callbacks\n",
        ")\n",
        "\n",
        "model = keras.models.load_model(\"jena_lstm.keras\")\n",
        "print(f\"Test MAE: {model.evaluate(test_dataset)[1]:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M6k-N5cDghiW"
      },
      "source": [
        "Plot learning curves:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pt-wCnf0giXo"
      },
      "outputs": [],
      "source": [
        "plot_loss_curve(history)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5IJbx_JRglCG"
      },
      "source": [
        "Maybe we need to regularise it. We add `recurrent_dropout`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "40fzRNfHgnZE"
      },
      "outputs": [],
      "source": [
        "inputs = keras.Input(shape=(sequence_length, raw_data.shape[-1]))\n",
        "x = keras.layers.LSTM(32, recurrent_dropout=0.25)(inputs)\n",
        "x = keras.layers.Dropout(0.5)(x)\n",
        "outputs = keras.layers.Dense(1)(x)\n",
        "model = keras.Model(inputs, outputs)\n",
        "callbacks = [\n",
        "    keras.callbacks.ModelCheckpoint(\"jena_lstm_dropout.keras\",\n",
        "    save_best_only=True)\n",
        "]\n",
        "\n",
        "model.compile(optimizer=\"rmsprop\", loss=\"mse\", metrics=[\"mae\"])\n",
        "history = model.fit(\n",
        "    train_dataset,\n",
        "    epochs=50,\n",
        "    validation_data=val_dataset,\n",
        "    callbacks=callbacks\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3amf0yyEg8p_"
      },
      "source": [
        "# Stacked LSTM"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2CpXSYv2g-qQ"
      },
      "outputs": [],
      "source": [
        "inputs = keras.Input(shape=(sequence_length, raw_data.shape[-1]))\n",
        "x = keras.layers.GRU(32, recurrent_dropout=0.5, return_sequences=True)(inputs)\n",
        "x = keras.layers.GRU(32, recurrent_dropout=0.5)(x)\n",
        "x = keras.layers.Dropout(0.5)(x)\n",
        "outputs = keras.layers.Dense(1)(x)\n",
        "model = keras.Model(inputs, outputs)\n",
        "callbacks = [\n",
        "    keras.callbacks.ModelCheckpoint(\"jena_stacked_gru_dropout.keras\",\n",
        "    save_best_only=True)\n",
        "]\n",
        "\n",
        "inputs = keras.Input(shape=(sequence_length, num_features))\n",
        "x = layers.LSTM(32, recurrent_dropout=0.2, unroll=True)(inputs)\n",
        "\n",
        "model.compile(optimizer=\"rmsprop\", loss=\"mse\", metrics=[\"mae\"])\n",
        "history = model.fit(train_dataset,\n",
        "    epochs=50,\n",
        "    validation_data=val_dataset,\n",
        "    callbacks=callbacks)\n",
        "model = keras.models.load_model(\"jena_stacked_gru_dropout.keras\")\n",
        "print(f\"Test MAE: {model.evaluate(test_dataset)[1]:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yAj6LZ76hVIt"
      },
      "source": [
        "Results:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HvJUbTUEhWUk"
      },
      "outputs": [],
      "source": [
        "plot_loss_curve(history)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9swpqAiJhYCf"
      },
      "source": [
        "# Try a bidirectional LSTM\n",
        "\n",
        "This one processes the input sequences twice: Once in chronological order, and once in reverse order.\n",
        "\n",
        "Does it help in predicting our weather data?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KKoz0LV5hc-Z"
      },
      "outputs": [],
      "source": [
        "inputs = keras.Input(shape=(sequence_length, raw_data.shape[-1]))\n",
        "x = keras.layers.Bidirectional(keras.layers.LSTM(16))(inputs)\n",
        "outputs = keras.layers.Dense(1)(x)\n",
        "model = keras.Model(inputs, outputs)\n",
        "\n",
        "\n",
        "model.compile(optimizer=\"rmsprop\", loss=\"mse\", metrics=[\"mae\"])\n",
        "history = model.fit(train_dataset,\n",
        "epochs=10,\n",
        "validation_data=val_dataset)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3t-IdnBIhpco"
      },
      "source": [
        "### <span style=\"color: red; font-weight: bold;\">Open exercises:<span>\n",
        "\n",
        "- Adjust the number of units in each recurrent layer in the stacked setup, as well as the amount of dropout. (The current choices are largely arbitrary and\n",
        "probably suboptimal.)\n",
        "\n",
        "- Adjust the learning rate used by the RMSprop optimizer, or try a different\n",
        "optimizer.\n",
        "\n",
        "- Try using a stack of Dense layers as the regressor on top of the recurrent layer, instead of a single Dense layer.\n",
        "\n",
        "- Improve the input to the model: try using longer or shorter sequences or a different sampling rate,"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "gpuType": "T4",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}