Copy of Surface_Observation_Prediction_in_Pytorch (1).ipynb

{"cells":[{"cell_type":"markdown","id":"df1dc0c0","metadata":{"id":"df1dc0c0"},"source":["# Surface Observation Postprocessing with Machine Learning\n","\n","[![](https://img.shields.io/badge/open-colab-yellow)](https://colab.research.google.com/github/ecmwf-projects/mooc-machine-learning-weather-climate/blob/main/tier_2/deep_learning/Surface_Observation_Prediction_in_Pytorch.ipynb) [![Gradient](https://assets.paperspace.io/img/gradient-badge.svg)](https://console.paperspace.com/github/ecmwf-projects/mooc-machine-learning-weather-climate/blob/main/tier_2/deep_learning/Surface_Observation_Prediction_in_Pytorch.ipynb) [![Open%20In%20SageMaker%20Studio%20Lab](https://studiolab.sagemaker.aws/studiolab.svg)](https://studiolab.sagemaker.aws/import/github/ecmwf-projects/mooc-machine-learning-weather-climate/blob/main/tier_2/deep_learning/Surface_Observation_Prediction_in_Pytorch.ipynb) \n","\n","Surface observation processing refers to the process of collecting and analyzing weather observations taken at the Earth's surface by weather stations and other observing platforms. This typically includes measuring temperature, humidity, wind speed and direction, precipitation, and other meteorological variables. The observations are then used to produce weather forecasts and warnings, as well as to analyze and understand weather patterns. The process typically involves quality control checks to ensure the observations are accurate and reliable, as well as the use of automated systems to process and disseminate the observations in near real-time.\n","\n","Surface observation processing is the process of collecting and analyzing weather information taken at ground level by weather stations and other observing platforms. These measurements include temperature, humidity, wind speed and direction, and precipitation. From this information, we can create weather forecasts and warnings and understand weather patterns.\n","\n","Machine learning can be applied in surface observation processing in several ways. One way is to use machine learning algorithms to help improve the accuracy of forecasts by analyzing large amounts of historical weather data. The algorithms can make more accurate predictions about future weather by learning from past weather patterns. Another way machine learning can be applied is in the quality control of the observations. Machine learning models can be trained to identify and flag any observations that may be incorrect or unreliable. This quality control can help ensure that only the most accurate and reliable weather information is used for forecasts and other purposes.\n","\n","In addition, Machine learning can be used to automate the process of surface observation processing, allowing the system to run more efficiently and quickly. This can help produce forecasts and warnings more quickly and with higher accuracy, providing more detailed and accurate information to the public.\n","\n","Overall, the use of Machine learning in surface observation processing can help improve the accuracy and speed of weather forecasting, and provide more detailed and accurate information to the public.\n","\n","In this notebooks, we are experimenting with deep learning techniques to solve a real-world problem that is relatively simple. The goal is to improve the temperature forecast of the ECMWF Integrated Forecast System (IFS) by using a technique called postprocessing.\n","\n","Specifically, we are trying to predict the difference between the temperature observed at a weather station and the temperature forecasted by the IFS at the closest point on the grid to the station. To achieve this, we will look at several different predictors based on physical principles.\n","\n","We will also get to explore several physically motivated predictors.\n","\n","_An example of using machine learning to model weather forecast error based on work of Matthew Chantry and Fenwick Cooper as part of work funded by [IFAB](https://www.ifabfoundation.org/). Modified by Jesper Dramsch for the [MOOC Machine Learning in Weather & Climate](https://ecmwf.int/mlwc-mooc) brough to you by ECMWF in partnership with IFAB._"]},{"cell_type":"markdown","id":"1dbf6081","metadata":{"id":"1dbf6081"},"source":["## Setting up the Notebook\n","\n","First, we need to install the necessary packages using `pip`.\n","\n","We are using Climetlab to load our prepared dataset, which we need to install\n","\n","### Data and Machine Learning Libraries\n","\n","Depending on the compute platform you are using, you may want to install the following packages in addition. These come pre-installed on some online environments like Google Colab, but not on all:\n","\n","- numpy\n","- matplotlib\n","- scikit-learn\n","- pytorch\n","- pytorch-lightning\n","\n","We added those to the `requirements.txt` for your convenience. If you're trying this exercise locally on Windows, it may be easier installing the libraries with `conda`. Pytorch has a guide to [install locally here](https://pytorch.org/get-started/locally/) and Tensorflow has a guide to [install locally here](https://www.tensorflow.org/install/pip), which may differ from they minimal code below."]},{"cell_type":"code","execution_count":null,"id":"d7439fdc","metadata":{"id":"d7439fdc"},"outputs":[],"source":["!pip install \"climetlab-mltc-surface-observation-postprocessing>=0.3.0\" --quiet\n","\n","# !pip install numpy matplotlib scikit-learn  --quiet\n","\n","# Pytorch\n","!pip install torch pytorch-lightning --quiet"]},{"cell_type":"markdown","id":"0724002c","metadata":{"id":"0724002c"},"source":["We will be using some tools that are commonly used in machine learning. \n","\n","Even though some of the mathematical concepts behind these tools are simple, it's often better to use pre-existing tools that have been tested and refined by others. This will save us time and make our analysis more reliable.\n","\n","<div class=\"alert alert-block alert-info\">\n","<b>Overview:</b> Numpy is the numerical Python library that is the de facto standard in dealing with scientific data in n-dimensional arrays\n","Matplotlib is a plotting library that can be used to create static, animated, and interactive visualizations of data.\n","Scikit-learn is a machine learning library for Python that provides simple and efficient tools for classic machine learning.\n","Tensorflow is a library for developing and deploying machine learning models, particularly neural networks.\n","Pytorch is a library for developing and deploying machine learning models, particularly deep learning models.</div>"]},{"cell_type":"code","execution_count":null,"id":"9716a266","metadata":{"id":"9716a266"},"outputs":[],"source":["import numpy as np\n","from sklearn import metrics\n","from sklearn.model_selection import train_test_split\n","from sklearn.preprocessing import StandardScaler\n","import matplotlib.pyplot as plt\n","\n","# CliMetLab will be used to manage the data loading.\n","import climetlab as cml\n","\n","# Import a tool for helping make figures\n","from climetlab_mltc_surface_observation_postprocessing.utils import imgBufferFromVectors\n"]},{"cell_type":"markdown","id":"c42b3279","metadata":{"id":"c42b3279"},"source":["## Loading the Data\n","\n","We will be loading a list that contains the errors of temperature forecasts made by the high-resolution forecast system of the European Centre for Medium-Range Weather Forecasts (ECMWF) for a 36-hour period. These forecasts will be compared with temperature measurements taken from about 8000 weather stations located around the world.\n","\n","<div class=\"alert alert-block alert-info\"><strong>Note:</strong> It's important to note that the dataset we're using has already gone through a crucial step of preparation known as data preprocessing. This means that any data that could be considered invalid or unreliable has already been removed. For example, data from weather stations with inconsistent measurement locations, duplicate values, or values that are physically impossible (such as temperatures over 100 degrees Celsius) have been removed. If you're starting a new project, this step of cleaning the data is essential.</div>"]},{"cell_type":"code","execution_count":null,"id":"a4d8e4da","metadata":{"id":"a4d8e4da"},"outputs":[],"source":["forecast_error = cml.load_dataset(\"mltc-surface-observation-postprocessing\", field=\"forecast_error\").to_numpy()\n"]},{"cell_type":"markdown","id":"28bbad8e","metadata":{"id":"28bbad8e"},"source":["We have a hypothesis that the errors in the forecast can be easily explained by the time of day the measurements were taken. To test this, we will load the local time of day for when the measurements were recorded."]},{"cell_type":"code","execution_count":null,"id":"d0b49363","metadata":{"id":"d0b49363"},"outputs":[],"source":["time_of_day = cml.load_dataset(\"mltc-surface-observation-postprocessing\", field=\"time_of_day\").to_numpy()\n"]},{"cell_type":"markdown","id":"e2ee1569","metadata":{"id":"e2ee1569"},"source":["Later we will add in the soil temperature, so we'll load this here to have the data loading step complete."]},{"cell_type":"code","execution_count":null,"id":"cb59e827","metadata":{"id":"cb59e827"},"outputs":[],"source":["soil_temperature = cml.load_dataset(\"mltc-surface-observation-postprocessing\", field=\"soil_temperature\").to_numpy()\n"]},{"cell_type":"markdown","id":"8418195e","metadata":{"id":"8418195e"},"source":["## Data Splitting\n","\n","We will be dividing the data used to make predictions and the data used to predict into two separate groups: training and testing. The training group will consist of 80% of the data and the testing group will consist of 20% of the data. This is done randomly, so that the data in each group is representative of the overall dataset.\n","\n","It's worth noting that the data has already been anonymized, so we can't use information about the weather stations or the time the measurements were taken to make sure that the data subsets are independent. If we had more metadata, we could use the methods we learned earlier to ensure that the data subsets are not dependent on each other.\n","\n","We use the function \"train_test_split\" from the scikit-learn library.\n","\n","The function takes an arbitrary number of inputs to split. Here we provide 3 inputs:\n","\n","- `forecast_error`: it is the data containing errors in forecasted temperatures,\n","- `time_of_day`: it is the data containing the local time of day when the measurements were taken,\n","- `soil_temperature`: it is the data containing the soil temperature at the location where the measurements were taken.\n","\n","The function also takes in two parameters:\n","\n","- `test_size`: it specifies what proportion of the data should be in the testing group. In this case, it is set to 0.2, which means that 20% of the data will be in the testing group and 80% of the data will be in the training group.\n","- `random_state`: it is used for initializing the random number generator and get the same split each time the script is run. Here it is set to `42`, which is an arbitrary number.\n","\n","The function returns 6 outputs:\n","\n","- `forecast_error_train`: it is the data containing errors in forecasted temperatures that will be used for training,\n","- `forecast_error_test`: it is the data containing errors in forecasted temperatures that will be used for testing,\n","- `time_of_day_train`: it is the data containing the local time of day when the measurements were taken that will be used for training,\n","- `time_of_day_test`: it is the data containing the local time of day when the measurements were taken that will be used for testing,\n","- `soil_temperature_train`: it is the data containing the soil temperature at the location where the measurements were taken that will be used for training,\n","- `soil_temperature_test`: it is the data containing the soil temperature at the location where the measurements were taken that will be used for testing.\n","\n","In summary, this code is splitting up the data into training and testing groups, using 80% of the data for training and 20% for testing. And it will be done in a random way and keep the same split each time the script is run by setting random_state to 42."]},{"cell_type":"code","execution_count":null,"id":"df69b984","metadata":{"id":"df69b984"},"outputs":[],"source":["(\n","    forecast_error_train,\n","    forecast_error_test,\n","    time_of_day_train,\n","    time_of_day_test,\n","    soil_temperature_train,\n","    soil_temperature_test,\n",") = train_test_split(\n","    forecast_error,\n","    time_of_day,\n","    soil_temperature,\n","    test_size=0.2,\n","    random_state=42,\n",")\n"]},{"cell_type":"markdown","id":"522f791e","metadata":{"id":"522f791e"},"source":["## Data Preprocessing\n","\n","The next step is to standardize the data. This is an important step that helps to prepare the data for training a neural network, especially in larger problems. To standardize the data, we will be using a method called `StandardScaler` method from scikit-learn library. This method divides each data point by the mean and standard deviation of the training dataset. By doing this, it ensures that the data has a mean of 0 and a standard deviation of 1.\n","\n","The code creates an instance of the `StandardScaler` class, and assigns it to the variable \"scaler\". Then we use the `fit()` method of this scaler object, which calculates the mean and variance of the `time_of_day_train` data. Note that we fit this object on the training data only, to avoid cheating and accidentally creating a model that seems like it works better than it actually does.\n","\n","We apply the transformation to `time_of_day_train` data and creates a new variable `time_of_day_train_norm` which contains the standardized data.\n","\n","The next line of code applies the same transformation to `time_of_day_test` data and creates a new variable `time_of_day_test_norm` which contains the standardized test data. Note that this only applies the standardization and doesn't learn from the test data, which is valid. This is one of the big benefits of using the `StandardScaler` object.\n","\n","In summary, this code is standardizing the `time_of_day` data by calculating the mean and variance of the training data and then dividing the training and testing data by those values resulting in a new set of data with mean = 0 and variance = 1.\n","\n","It's worth noting that there are other ways to prepare data for a neural network, and sometimes a different approach may be more appropriate. For example, when working with precipitation data, it may be more appropriate to use a log transform instead of standardizing the data."]},{"cell_type":"code","execution_count":null,"id":"b3ffb07a","metadata":{"id":"b3ffb07a"},"outputs":[],"source":["scaler = StandardScaler()\n","scaler.fit(time_of_day_train)\n","print(\"Mean:\", scaler.mean_, \", Variance:\", scaler.var_)\n","\n","time_of_day_train_norm = scaler.transform(time_of_day_train)\n","time_of_day_test_norm = scaler.transform(time_of_day_test)\n"]},{"cell_type":"markdown","id":"58759cc8","metadata":{"id":"58759cc8"},"source":["### Pytorch DataLoaders\n","\n","Pytorch expects a Dataloader object to provide access to the data we just created.\n","\n","There are many Dataset and Dataloader objects pre-defined in the Pytorch ecosystem and we'll use the `TensorDataset` to wrap our data to pass into a `DataLoader`.\n","\n","But first we import the necessary parts from the Pytorch library:"]},{"cell_type":"code","execution_count":null,"id":"68672bac","metadata":{"id":"68672bac"},"outputs":[],"source":["import torch\n","import pytorch_lightning as pl\n","\n","from torch import nn\n","from torch.nn import functional as F\n","from torch.optim import Adam\n","from torch.utils.data import TensorDataset, DataLoader\n"]},{"cell_type":"markdown","id":"e53deb0b","metadata":{"id":"e53deb0b"},"source":["The first step is to convert the NumPy arrays to PyTorch tensors. This is done using the `torch.from_numpy()` function, which converts a NumPy array to a PyTorch tensor. The `.float()` at the end of the line is used to convert the tensors to a float data type.\n","\n","The next step is to create a `TensorDataset` from the data. This dataset is created by passing the data and labels tensors to the TensorDataset constructor.\n","\n","Finally, we create a `DataLoader` from the dataset. A DataLoader is a PyTorch utility that allows us to iterate over the dataset in batches. It takes the dataset and some options as input, such as the batch size and whether to shuffle the data. In this case, the batch size is set to 32, and the data is set to be shuffled before each epoch.\n","\n","With this, you have a DataLoader that can be used to iterate over the data in batches during the training process. The DataLoader is useful for loading large datasets that do not fit into memory all at once, as it loads only a batch at a time.\n","\n","## Validation Strategy\n","\n","We learned that we can only use the test data for a final evaluation of the model performance after we are done tweaking the model.\n","\n","Therefore we need to split the training data again and define a validation dataloader as well."]},{"cell_type":"code","execution_count":null,"id":"488c275f","metadata":{"id":"488c275f"},"outputs":[],"source":["class SOPDataModule(pl.LightningDataModule):\n","    def __init__(self, train_data, train_labels, test_data, test_labels, random_state=None):\n","        super().__init__()\n","        \n","        self.train_data, self.val_data, self.train_label, self.val_label = train_test_split(\n","           train_data, train_labels, test_size=0.15, random_state=random_state\n","        )\n","        self.test_data = test_data\n","        self.test_label = test_labels\n","\n","        self.batch_size = 128\n","    \n","    def prepare_data(self):\n","        self.dataset_train = TensorDataset(torch.from_numpy(self.train_data).float(), torch.from_numpy(self.train_label).float())\n","        self.dataset_val = TensorDataset(torch.from_numpy(self.val_data).float(), torch.from_numpy(self.val_label).float())\n","        self.dataset_test = TensorDataset(torch.from_numpy(self.test_data).float(), torch.from_numpy(self.test_label).float())\n","\n","    def train_dataloader(self):\n","        return DataLoader(self.dataset_train, batch_size=self.batch_size)\n","\n","    def val_dataloader(self):\n","        return DataLoader(self.dataset_val, batch_size=self.batch_size)\n","\n","    def test_dataloader(self):\n","        return DataLoader(self.dataset_test, batch_size=self.batch_size)\n"]},{"cell_type":"code","execution_count":null,"id":"31aa8aa9","metadata":{"id":"31aa8aa9"},"outputs":[],"source":["datamodule = SOPDataModule(time_of_day_train_norm, forecast_error_train, time_of_day_test_norm, forecast_error_test)"]},{"cell_type":"markdown","id":"1cd5a3ea","metadata":{"id":"1cd5a3ea"},"source":["## Model Building\n","\n","In Pytorch Lightning we commonly build models as a class.\n","\n","Within these classes, we can define models in multiple ways\n","\n","This code defines two functions: \"fullyconnected_sequential\" and \"fullyconnected_functional\". Both functions create a model with a similar structure: it starts with an input layer, followed by several dense layers, and ends with an output layer. \n"]},{"cell_type":"markdown","id":"a3c7f4f4","metadata":{"id":"a3c7f4f4"},"source":["### Sequential Model\n","\n","This code is defining a PyTorch Lightning module called `FullyConnected`. PyTorch Lightning is a high-level library for PyTorch that makes it easier to train and evaluate deep learning models by providing a simple and consistent interface.\n","\n","The constructor takes several parameters:\n","\n","- `input_shape`: an integer representing the number of predictors.\n","- `width`: an integer representing how wide the layers should be.\n","- `depth`: an integer representing how many layers the network should have.\n","- `activation`: a string representing the non-linear activation to use (e.g. 'relu', 'sigmoid', etc.).\n","- `learning_rate`: a float representing the learning rate of the optimizer (default value is 10^-3).\n","- `final_activation`: a string representing the activation function for the final layer(default value is 'linear').\n","\n","We create a number of consecutive linear layers with the activation and width as input and output size respectively. \n","\n","`nn.Sequential` is used as the model, which is a container that sequentially applies a number of modules.\n","\n","Then it adds the final activation and a linear layer as the output layer with the width and 1 as input and output size respectively.\n","\n","The `training_step()` function is used to define the training loop. It takes a batch of data and the batch index, then it passes the data through the model, calculates the loss using the mean squared error loss function `F.mse_loss`, and logs the loss to the training log.\n","\n","The `configure_optimizers` function is used to define the optimizer used to train the model. In this case, it is using the `Adam` optimizer with the learning rate that is passed to the class constructor.\n","\n","This code provides a basic example of how to create a PyTorch Lightning module, define the model architecture and the training loop, and how to configure the optimizer."]},{"cell_type":"code","execution_count":null,"id":"f7894d4d","metadata":{"id":"f7894d4d"},"outputs":[],"source":["class FullyConnectedSequential(pl.LightningModule):\n","    def __init__(self, input_shape, width, depth, activation=None, learning_rate=10**(-3), final_activation=None):\n","        super().__init__()\n","        \n","        self.input_shape = input_shape\n","        self.width = width\n","        self.depth = depth\n","\n","        self.learning_rate = learning_rate\n","\n","        self.activation = activation() if activation is not None else nn.ReLU()\n","        self.final_activation = final_activation() if final_activation is not None else nn.Identity()\n","\n","        self.hidden = []\n","        for _ in range(depth):\n","            self.hidden.append(self.activation)\n","            self.hidden.append(nn.Linear(width, width))\n","\n","        self.model = nn.Sequential(nn.Linear(input_shape, width), *self.hidden, self.final_activation, nn.Linear(width, 1))\n","\n","    def forward(self, x):\n","        return self.model(x)\n","\n","    def training_step(self, batch, batch_idx):\n","        x, y = batch\n","        y_hat = self(x)\n","        loss = F.mse_loss(y_hat, y)\n","        self.log(\"train_loss\", loss)\n","        return loss\n","    \n","    def validation_step(self, batch, batch_idx):\n","        x, y = batch\n","        y_hat = self(x)\n","        loss = F.mse_loss(y_hat, y)\n","        self.log(\"validation_loss\", loss)\n","        return loss\n","\n","    def test_step(self, batch, batch_idx):\n","        x, y = batch\n","        y_hat = self(x)\n","        loss = F.mse_loss(y_hat, y)\n","        self.log(\"test_loss\", loss)\n","        return loss\n","\n","    def configure_optimizers(self):\n","        opt = Adam(self.parameters(), lr=self.learning_rate)\n","        return opt\n"]},{"cell_type":"markdown","id":"a5b74335","metadata":{"id":"a5b74335"},"source":["### Functional Model\n","\n","This code is defining a PyTorch Lightning module called FullyConnectedFunctional. This module is similar to the FullyConnected module that you previously explained, but it uses a functional approach to define the model architecture.\n","\n","The constructor takes several parameters:\n","\n","- `input_shape`: an integer representing the number of predictors.\n","- `width`: an integer representing how wide the layers should be.\n","- `depth`: an integer representing how many layers the network should have.\n","- `activation`: a string representing the non-linear activation to use (e.g. 'relu', 'sigmoid', etc.).\n","- `learning_rate`: a float representing the learning rate of the optimizer (default value is 10^-3).\n","- `final_activation`: a string representing the activation function for the final layer(default value is 'linear').\n","\n","The input, hidden and output layers are defined as functions, each function takes in a tensor and applies the linear layer and activation functions to it and return the output tensor.\n","\n","The `forward` function is used to define the forward pass of the model. It receives an input tensor and applies the input layer to it. Then it iteratively applies the hidden layer a number of times equal to the depth of the model. Finally, it applies the output layer and return the output tensor.\n","\n","The `training_step` function is used to define the training loop. It takes a batch of data and the batch index, then it passes the data through the model, calculates the loss using the mean squared error loss function (F.mse_loss), and logs the loss to the training log.\n","\n","The `configure_optimizers` function is used to define the optimizer used to train the model. In this case, it is using the `Adam` optimizer with the learning rate that is passed to the class constructor.\n","\n","This code provides a basic example of how to create a PyTorch Lightning module using a functional approach, define the model architecture and the training loop, and how to configure the optimizer."]},{"cell_type":"code","execution_count":null,"id":"848ff69b","metadata":{"id":"848ff69b"},"outputs":[],"source":["\n","class FullyConnectedFunctional(pl.LightningModule):\n","    def __init__(self, input_shape, width, depth, activation=None, learning_rate=10**(-3), final_activation=None):\n","        super().__init__()\n","        \n","        self.input_shape = input_shape\n","        self.width = width\n","        self.depth = depth\n","        \n","        self.learning_rate = learning_rate\n","\n","        self.activation = activation() if activation is not None else nn.ReLU()\n","        self.final_activation = final_activation() if final_activation is not None else nn.Identity()\n","        \n","        self.input = nn.Linear(input_shape, width)\n","        \n","        for i in range(self.depth):\n","            setattr(self, f\"hidden{i}\", nn.Linear(width, width))\n","        self.output = nn.Linear(width, 1)\n","\n","    def forward(self, x):\n","        x = self.input(x)\n","        for i in range(self.depth):\n","            layer = getattr(self, f\"hidden{i}\")\n","            x = layer(self.activation(x))\n","        x = self.final_activation(x)\n","        return self.output(x)\n","\n","    def training_step(self, batch, batch_idx):\n","        x, y = batch\n","        y_hat = self(x)\n","        loss = F.mse_loss(y_hat, y)\n","        self.log(\"train_loss\", loss)\n","        return loss\n","    \n","    def validation_step(self, batch, batch_idx):\n","        x, y = batch\n","        y_hat = self(x)\n","        loss = F.mse_loss(y_hat, y)\n","        self.log(\"validation_loss\", loss)\n","        return loss\n","    \n","    def test_step(self, batch, batch_idx):\n","        x, y = batch\n","        y_hat = self(x)\n","        loss = F.mse_loss(y_hat, y)\n","        self.log(\"test_loss\", loss)\n","        return loss\n","\n","    def configure_optimizers(self):\n","        opt = Adam(self.parameters(), lr=self.learning_rate)\n","        return opt\n"]},{"cell_type":"markdown","id":"14f3d60b","metadata":{"id":"14f3d60b"},"source":["#### Adam Optimiser\n","\n","When training a neural network, we use an optimization algorithm to adjust the parameters of the model in order to minimize the loss function. Notice that the optimiser we use is called `Adam` instead of the classic gradient descent algorithm described during the lessons.\n","\n","Gradient descent is a simple optimization algorithm that updates the model's parameters by taking a step in the direction of the negative gradient of the loss function. The step size is determined by a learning rate parameter, which controls how large the step should be. Gradient descent can work well for simple models, but can be slow to converge for more complex models.\n","\n","Adam, on the other hand, is a more advanced optimization algorithm that is based on gradient descent. It uses an adaptive learning rate, which means that the step size is adjusted for each parameter based on the historical gradient information. This can help the optimization converge faster and more efficiently. Adam also includes additional features such as momentum and RMSprop, which can further improve the optimization process. This makes Adam more similar to the Conjuage Gradient method in classic numerical optimisation. Choosing Adam can improve the optimization process and make it converge faster and more efficiently.\n","\n","### Creating the Model\n","Create neural network object you can use either of the functions we defined above\n","\n","\n","Try changing the depth and see what the summary says:\n","\n","<div class=\"alert alert-block alert-warning\"> \n","<b>Alternative:</b><code>\n","model = FullyConnectedFunctional(\n","    time_of_day_train_norm.shape[1], depth=3, width=32, activation=\"tanh\", learning_rate=10**(-3)\n",")</code>\n","</div>"]},{"cell_type":"code","execution_count":null,"id":"b717e3fd","metadata":{"id":"b717e3fd"},"outputs":[],"source":["model = FullyConnectedSequential(\n","    time_of_day_train_norm.shape[1], depth=3, width=32, activation=nn.Tanh, learning_rate=10 ** (-3)\n",")\n","\n","pl.utilities.model_summary.ModelSummary(model, max_depth=-1)\n"]},{"cell_type":"markdown","id":"34989018","metadata":{"id":"34989018"},"source":["<div class=\"alert alert-block alert-info\">\n","<b>Note:</b> We have 2000 parameters to learn. This is many fewer than the number of examples, so we should not be too worried about overfitting.</div>\n"]},{"cell_type":"markdown","id":"4d17e383","metadata":{"id":"4d17e383"},"source":["## Model training\n","\n","The `fit` function implements the training loop for a model in Pytorch lightning. \n","\n","The first two arguments passed to the `fit` function is the model. In this case, the training data is `time_of_day_train_norm` and the labels are `forecast_error_train`. These are the inputs and outputs that the model will use to learn during training as supplied by the dataloader.\n","\n","The `epochs` argument is the number of times the model will cycle through the data.\n","\n","Once the training is done, the model will be optimized to make predictions on new data with the same features."]},{"cell_type":"code","execution_count":null,"id":"40b46cd2","metadata":{"id":"40b46cd2"},"outputs":[],"source":["try:\n","    trainer = pl.Trainer(max_epochs=3, accelerator='gpu', devices=1)\n","except:\n","    print(\"Running on CPU will be significantly slower. Try activating a GPU backend.\")\n","    trainer = pl.Trainer(max_epochs=3)\n","\n","trainer.fit(model, datamodule)\n"]},{"cell_type":"markdown","id":"a7f06c96","metadata":{"id":"a7f06c96"},"source":["In other problems it may take far longer training to reach this point, or even be impossible if a network is not appropriate for the data provided.\n","\n","## Model Evaluation\n","\n","Finally we can do a numerical evaluation how well the model performs on unseen data.\n","\n","Purely, as a teaching device, we'll also evaluate the model on training data. This number is relatively meaningless standing by itself, as the model will always perform on training data as long as the model converged. The validation data is unseen as training, but was involved in the optimization process, which can have indirect optimization consequences, by us adjusting the architecture manually and generally making design choices based on the validation results.\n","\n","But it's good to learn, why we test on actual unseen data and use the training value as the lower bound for our expectations of the final evaluation of the test data."]},{"cell_type":"code","execution_count":null,"id":"cbbabb13","metadata":{"id":"cbbabb13"},"outputs":[],"source":["train_loss = trainer.test(dataloaders=datamodule.train_dataloader())\n","print(\"Train loss:\", train_loss)\n","\n","val_loss = trainer.test(dataloaders=datamodule.val_dataloader())\n","print(\"Validation loss:\", val_loss)\n"]},{"cell_type":"markdown","id":"583c3d36","metadata":{"id":"583c3d36"},"source":["Here we see the actual performance of our model compared based on true unseen data."]},{"cell_type":"code","execution_count":null,"id":"18b86546","metadata":{"id":"18b86546"},"outputs":[],"source":["test_loss = trainer.test(dataloaders=datamodule.test_dataloader())\n","print(\"Test loss:\", test_loss)\n"]},{"cell_type":"markdown","id":"7091771c","metadata":{"id":"7091771c"},"source":["### Visualization for Evaluation\n","\n","We can create some fake data to cover the entire day and get predictions to visualize the entire space of the data. \n","\n","This is not always possible, but that's what we have this tutorial data for to gain an intuitive understanding of the methods we apply. This creates an array of time of day values spanning from 0 to 24, with a step of 0.001. This is a full range of time of day values that we want to make predictions on. The array is then reshaped to have a single axis, this is done with the `[..., np.newaxis]` slice.\n","\n","Then we apply the same standardization procedure (using the fitted `scaler` object) on the full range of time of day values as the one used on the training and test data, to ensure that the new data has the same scale as the data the model was trained on.\n","\n","Finally, the code uses the in-built model's function to make predictions on the normalized full range of time of day values. The predict function returns an array of predicted forecast errors, one for each input time of day value. The variable `forecast_error_full_range_predicted` stores these prediction."]},{"cell_type":"code","execution_count":null,"id":"cd3ae3bf","metadata":{"id":"cd3ae3bf"},"outputs":[],"source":["time_of_day_full_range = np.arange(0, 24, 0.001)[..., np.newaxis]\n","time_of_day_full_range_norm = scaler.transform(time_of_day_full_range)\n","forecast_error_full_range_predicted = model(torch.Tensor(time_of_day_full_range_norm)).detach().numpy()\n"]},{"cell_type":"markdown","id":"ddad7ab6","metadata":{"id":"ddad7ab6"},"source":["Then we create an image of time of day values against forecast error values. The image is created by using a function called `imgBufferFromVectors` which takes a set of `x` and `y` values and creates an image.\n","\n","The function takes these arguments:\n","\n","- `time_of_day_test`: the x values, representing the time of day\n","- `forecast_error_test`: the y values, representing the forecast error\n","- `nx`: the number of pixels in the x-axis of the image\n","- `ny`: the number of pixels in the y-axis of the image\n","- `extent`: the range of values to be plotted in the image, if empty it will use the min and max of the x and y values.\n","- `calc_average`: a boolean flag indicating whether to calculate the average of the input values or not.\n","\n","The function then returns the image as a buffer, the extent of the axis and the count of values in the image."]},{"cell_type":"code","execution_count":null,"id":"7987db66","metadata":{"id":"7987db66"},"outputs":[],"source":["tod_buffer, ax_extent, count = imgBufferFromVectors(\n","    time_of_day_test, forecast_error_test, nx=256, ny=256, extent=[], calc_average=False\n",")\n"]},{"cell_type":"markdown","id":"d657e2f7","metadata":{"id":"d657e2f7"},"source":["Let's plot the number of measurements at each time of day against the forecast error.\n","\n","The red line shows the best fit.\n","\n","<div class=\"alert alert-block alert-info\"> <b>Note:</b> The model (red line) is better than nothing, but does not account for the variation in 2m temperature error at single time of day (eg. at 14:00). Another predictor is required. We try using the soil temperature. In comparison to the regression model, the result of the neural etwork is less smooth which can be a good OR a bad thing.</div>"]},{"cell_type":"code","execution_count":null,"id":"eade6b39","metadata":{"id":"eade6b39"},"outputs":[],"source":["plt.imshow(np.log((count == 0.0) + count), cmap=\"Blues\", origin=\"lower\", extent=ax_extent, aspect=\"auto\")\n","\n","plt.xlim([0, 24])\n","plt.grid()\n","plt.xlabel(\"Time of day\")\n","plt.ylabel(\"2m temperature error ($^\\mathrm{o}$C)\")\n","cb = plt.colorbar()\n","cb.set_label(\"Log( number of measurements )\")\n","\n","# Line of best fit\n","plt.plot(time_of_day_full_range, forecast_error_full_range_predicted, \"red\")\n","\n","plt.show()\n"]},{"cell_type":"markdown","id":"3d2794b5","metadata":{"id":"3d2794b5"},"source":["## Improving the Model\n","\n","Now we can start to play around.\n","\n","Let's add more predictors. \n","\n","We will add a second predictor, the model soil temperature. As we will see below, there is a complex structure in the pattern of forecast error in the 2D space of time of day & soil temperature. \n","\n","How well can the model learn this?\n","\n","### Visualize the Soil Temperature Data\n","\n","This code is plotting a 2D histogram of the relationship between forecast error, soil temperature and time of day. \n","\n","Notice how blue shows that the observation is warmer than the forecast and red shows that forecast is warmer than the observation."]},{"cell_type":"code","execution_count":null,"id":"b379f24d","metadata":{"id":"b379f24d"},"outputs":[],"source":["# Make image of the error with the new predictor\n","buffer, ax_extent, count = imgBufferFromVectors(\n","    soil_temperature_test, time_of_day_test, forecast_error_test, 128, 256, extent=[], calc_average=True\n",")\n","\n","# Plot the image of the error\n","plt.imshow(buffer, vmin=-5, vmax=5, cmap=\"seismic\", origin=\"lower\", extent=ax_extent, aspect=\"auto\")\n","\n","plt.grid()\n","plt.xlabel(\"Soil temperature ($^o$C)\")\n","plt.ylabel(\"Local time of day (hours)\")\n","cb = plt.colorbar()\n","cb.set_label(\"Forecast - Observation ($^o$C)\")\n","\n","plt.show()\n"]},{"cell_type":"markdown","id":"e2055228","metadata":{"id":"e2055228"},"source":["In the image above we see clear structure, where the in evening the forecast is too warm if the soil is frozen and too cold if the soil is not.\n","\n","Can we learn a good representation of this error pattern?\n","\n","### Training a new model\n","\n","We'll train a neural network model using two predictors, time of day and soil temperature. \n","\n","The data for the two predictors is first concatenated and then standardized using `StandardScaler` again. Then we construct a functional model with depth and width of 3 and 32 respectively, using `tanh` activation function to mix things up and train.\n","\n","<div class=\"alert alert-block alert-warning\">Again, this may take a couple of minutes</div>"]},{"cell_type":"code","execution_count":null,"id":"285729c8","metadata":{"id":"285729c8"},"outputs":[],"source":["# Create the input features\n","X_train = np.concatenate([time_of_day_train, soil_temperature_train], axis=-1)\n","X_test = np.concatenate([time_of_day_test, soil_temperature_test], axis=-1)\n","\n","# Standardise the data\n","scaler = StandardScaler()\n","scaler.fit(X_train)\n","X_train = scaler.transform(X_train)\n","X_test = scaler.transform(X_test)\n","\n","# Construct the Data Module\n","datamodule = SOPDataModule(X_train, forecast_error_train, X_test, forecast_error_test)\n","\n","# Construct a model, this time with 2 predictors\n","model_with_soil = FullyConnectedFunctional(\n","    X_train.shape[-1], depth=3, width=32, activation=nn.Tanh, learning_rate=10 ** (-3)\n",")\n","\n","# Alternative activation functions are tanh, sigmoid, softmax, relu, softplus...\n","try:\n","    trainer = pl.Trainer(max_epochs=3, accelerator='gpu', devices=1)\n","except:\n","    print(\"Training on CPU may be significantly slower. Try activating a GPU backend.\")\n","    trainer = pl.Trainer(max_epochs=3)\n","\n","trainer.fit(model_with_soil, datamodule)\n"]},{"cell_type":"markdown","id":"aa19cd5f","metadata":{"id":"aa19cd5f"},"source":["### Evaluating the New Model\n","\n","Let's do some evaluation on the unseen test data!\n","\n","We'll calculate the mean absolute error and the RMSE of the uncorrected forecast and the corrected forecast."]},{"cell_type":"code","execution_count":null,"id":"a0baf12e","metadata":{"id":"a0baf12e"},"outputs":[],"source":["# Calculate Root Mean Square error of predictions:\n","\n","zero_test = 0.0 * forecast_error_test\n","print(\"Mean Absolute Error Uncorrected:\", metrics.mean_absolute_error(zero_test, forecast_error_test))\n","print(\"Root Mean Squared Error Uncorrected:\", np.sqrt(metrics.mean_squared_error(zero_test, forecast_error_test)))\n","\n","forecast_corrected = forecast_error_test - model_with_soil(torch.Tensor(X_test)).detach().numpy()\n","\n","print(\"Mean Absolute Error Corrected:\", metrics.mean_absolute_error(zero_test, forecast_corrected))\n","print(\"Root Mean Squared Error Corrected:\", np.sqrt(metrics.mean_squared_error(zero_test, forecast_corrected)))\n"]},{"cell_type":"markdown","id":"af4c9c2a","metadata":{"id":"af4c9c2a"},"source":["By using two predictors we can impove the forecast of 2m-temperature by ~0.1C.\n","\n","#### Visual Evaluation\n","\n","Let's see how the model corrects the forecast across the full 2D space.\n","\n","We'll run the trained model over a range of values for the input variables `soil_temperature` and `time_of_day`. \n","\n","It first defines the range of values for these variables using `nx` and `ny`, which are the dimensions of the input buffer, and the `ax_extent` which is the range of values observed in the data. Then it creates the input data by stacking the meshgrid of these variable ranges using `np.meshgrid`, and transforms this input data using the `scaler` that was fit earlier on the training data. Then we predict the output using this input data and reshapes this output back to 2D plot using `raw_pred.reshape(input_buffer.shape[:-1])` for us to plot!"]},{"cell_type":"code","execution_count":null,"id":"68ffa23d","metadata":{"id":"68ffa23d"},"outputs":[],"source":["# Run the fit model over the plot domain\n","\n","# The x and y values of each point in the plot image, this covers\n","# the range of data observed in the data\n","nx = buffer.shape[0]\n","ny = buffer.shape[1]\n","x_st = np.linspace(ax_extent[0], ax_extent[1], nx)  # Represents soil_temperature\n","y_tod = np.linspace(ax_extent[2], ax_extent[3], ny)  # Represents time_of_day\n","\n","# Create the input data\n","input_buffer = np.stack(np.meshgrid(y_tod, x_st), axis=-1)\n","X_plot = scaler.transform(input_buffer.reshape(-1, X_train.shape[-1]))\n","\n","# Predict and reshape prediction back to 2D plot\n","raw_pred = model_with_soil(torch.Tensor(X_plot)).detach().numpy()\n","model_buffer = raw_pred.reshape(input_buffer.shape[:-1])\n"]},{"cell_type":"markdown","id":"953545c5","metadata":{"id":"953545c5"},"source":["Let's compare the images!\n","\n","Notice that we're fixing the range of the plot between -5 and 5 to make these plots intuitively comparable."]},{"cell_type":"code","execution_count":null,"id":"13f0ab18","metadata":{"id":"13f0ab18"},"outputs":[],"source":["# Plot the image of the error\n","plt.imshow(buffer, vmin=-5, vmax=5, cmap=\"seismic\", origin=\"lower\", extent=ax_extent, aspect=\"auto\")\n","\n","plt.grid()\n","plt.title(\"Forecast Error or Raw Data\")\n","plt.xlabel(\"Soil temperature ($^o$C)\")\n","plt.ylabel(\"Local time of day (hours)\")\n","cb = plt.colorbar()\n","cb.set_label(\"Forecast - Observation ($^o$C)\")\n","\n","plt.show()\n","\n","# Plot the model output where we have data\n","plt.imshow(\n","    (model_buffer) * (count > 0), vmin=-5, vmax=5, cmap=\"seismic\", origin=\"lower\", extent=ax_extent, aspect=\"auto\"\n",")\n","\n","plt.grid()\n","plt.title(\"Model Output of Predicted Forecast Error\")\n","plt.xlabel(\"Soil temperature ($^o$C)\")\n","plt.ylabel(\"Local time of day (hours)\")\n","cb = plt.colorbar()\n","cb.set_label(\"Neural Network model ($^o$C)\")\n","\n","plt.show()\n","\n","# Plot the model over the whole domain\n","plt.imshow(model_buffer, vmin=-5, vmax=5, cmap=\"seismic\", origin=\"lower\", extent=ax_extent, aspect=\"auto\")\n","\n","plt.grid()\n","plt.title(\"Model Output of Extrapolating to the Full Domain\")\n","plt.xlabel(\"Soil temperature ($^o$C)\")\n","plt.ylabel(\"Local time of day (hours)\")\n","cb = plt.colorbar()\n","cb.set_label(\"Neural Network model ($^o$C)\")\n","\n","plt.show()\n"]},{"cell_type":"markdown","id":"782f5fc6","metadata":{"id":"782f5fc6"},"source":["## Conclusion\n","\n","As we add more predictors and complexity to our model, it becomes better and better. But, our measurements of the forecast error are not perfect. There is some noise. We don't want our model to capture this. But an overly complex model will. We need a simpler model, or more data. Fitting a model that is too complex for the data is causes overfitting. \n","\n","Here we have a large number of observations and still relatively small number of free parameters, so overfitting is unlikely. We see that during training our errors on the training & testing dataset are comparable.\n","\n","Away from where data has been provided the model does not have constraints. We should not trust this part of the space without additional validation.\n","\n","We also see different model corrections for hour 23 and hour 0, when these should be closely correlated, as they're mere seconds away from each other.\n"]},{"cell_type":"markdown","id":"f4e7c3e0","metadata":{"id":"f4e7c3e0"},"source":["# Exercises\n","\n","1. Can you beat these predictions by changing the model architecture?\n","\n","Play around with activations, depth, width, learning rate of the neural network.\n","\n","2. Is there a way of building in any prior knowledge to even this simple setup?\n","\n","Perhaps you can encode the fact that 0 hour follows 24? Does this help the prediction?"]},{"cell_type":"markdown","id":"d833525e","metadata":{"id":"d833525e"},"source":["Now that you've had a proper introduction to deep learning, you may quench your curiosity and go back to the Weatherbench notebook from tier 1 to inspect a convolutional neural network.\n"]}],"metadata":{"accelerator":"GPU","colab":{"provenance":[{"file_id":"https://github.com/ecmwf-projects/mooc-machine-learning-weather-climate/blob/main/tier_2/deep_learning/Surface_Observation_Prediction_in_Pytorch.ipynb","timestamp":1681343881333}]},"gpuClass":"standard","kernelspec":{"display_name":"Python 3 (ipykernel)","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.9.13"},"vscode":{"interpreter":{"hash":"f49206fcf84a9145e7e21228cbafa911d1ac18292303b01e865d8267a9c448f7"}}},"nbformat":4,"nbformat_minor":5}