In [27]:
# Imports
# Note that for the theano model to work, the keras theano backend must be used!
from importlib import reload
import emos_network_theano; reload(emos_network_theano)
from emos_network_theano import EMOS_Network
from losses import crps_cost_function
import utils; reload(utils)
from utils import *
import keras_models; reload(keras_models)
from keras_models import *
In [2]:
# Basic setup
DATA_DIR = '/Volumes/STICK/data/ppnn_data/' # Mac
# DATA_DIR = '/project/meteo/w2w/C7/ppnn_data/' # LMU
results_dir = '../results/' # Where to save post-processed predictions
window_size = 25 # Days in rolling window
fclt = 48 # Forecast lead time in hours
The interpolated raw data set contains the observations and forecasts from the 50 members for each station and each day. Using the get_train_test_sets function we will convert this raw dataset to a format suitable for the networks.
Here we pick a forecast date and return as the training data all previous days within the window size previous to the start of the forecast. The test data is simply the data for the chosen forecast date. The 50 member ensemble is summarized by the mean and the standard deviation. Additionally, we remove all data where the observations are missing. Finally, the inputs are scaled. For this we simply divide each feature by its maximum in the training set.
In [3]:
# Chose a random date to illustrate the algorithm
date_str = '2011-02-14' # This is our standard date format
In [4]:
# Load training and test set
train_set, test_set = get_train_test_sets(DATA_DIR, predict_date=date_str,
fclt=fclt, window_size=window_size)
These two sets are objects which contain all the data and some meta information.
In [5]:
train_set.feature_names
Out[5]:
In [6]:
train_set.features.shape, train_set.targets.shape, train_set.date_strs.shape
Out[6]:
In [7]:
test_set.features.shape, test_set.targets.shape, test_set.date_strs.shape
Out[7]:
To start with we build the model in pure theano. The model is defined in a separate script EMOS_network_theano.py. The network uses a custom CRPS loss function which is defined in crps_loss.py.
The EMOS_Network class is build to work in a similar way to keras models. For the fitting we are using gradient descent. Since we are using the entire dataset for each update, it is not stochastic. An early stopping algorithm is built into the fitting function. It stops training if the average training CRPS of the last 5 steps is decreasing by less than a parameter delta.
In [8]:
# Define some model parameters
lr = np.asarray(0.1, dtype='float32') # The learning rate
early_stopping_delta = 1e-4 # How much the CRPS must improve before stopping
steps_max = 1000 # How many steps to fit at max
In [9]:
# Set up the theano model
model_theano = EMOS_Network()
In [10]:
# Split the features into means and standard deviation
train_set.feature_names
Out[10]:
In [11]:
train_mean = train_set.features[:, 0]
train_std = train_set.features[:, 1]
test_mean = test_set.features[:, 0]
test_std = test_set.features[:, 1]
In [12]:
# Train the model for some steps
model_theano.fit(train_mean, train_std, train_set.targets, steps_max,
(test_mean, test_std, test_set.targets), lr=lr,
early_stopping_delta=early_stopping_delta)
# Output is the training CRPS and the test CRPS
Out[12]:
To compare the network model with the standard EMOS we will run it from 1 January 2016 to 31 December 2016. When looping over the days we are not resetting the model weights for each day. This drastically reduces the training time with identical results.
In [13]:
# Get start and stop indices
date_str_start = '2016-01-01'
date_str_stop = '2017-01-01'
In [14]:
model_theano = EMOS_Network()
In [17]:
# This function loops over the days.
train_crps_list, valid_crps_list, results_df = loop_over_days(
DATA_DIR,
model_theano,
date_str_start, date_str_stop,
window_size=window_size,
fclt=fclt,
epochs_max=steps_max,
early_stopping_delta=early_stopping_delta,
lr=lr,
verbose=0,
model_type='EMOS_Network_theano')
Note that before restructuring the data preparation function, this was significantly slower. It might have to do with the creation of the meta data arrays.
In [18]:
# Let's see what the mean prediction CRPS is
np.mean(valid_crps_list)
Out[18]:
In [19]:
# Save the results
results_df.to_csv(results_dir + 'emos_network_rolling_window.csv')
This file is then read by the evaluation script.
In [20]:
model_keras = build_EMOS_network_keras(compile=True, optimizer='sgd', lr=0.1)
In [21]:
model_keras.summary()
In [22]:
# This way we have the gradient descent on the whole training set just as in theano
batch_size = train_mean.shape[0]
batch_size
Out[22]:
In [23]:
model_keras.fit([train_mean, train_std], train_set.targets, epochs=steps_max,
batch_size=batch_size,
validation_data=[[test_mean, test_std], test_set.targets],
verbose=0,
callbacks=[EarlyStopping(monitor='loss',
min_delta=early_stopping_delta,
patience=2)]);
Out[23]:
In [24]:
# Get train and test CRPS
(model_keras.evaluate([train_mean, train_std], train_set.targets, batch_size, verbose=0),
model_keras.evaluate([test_mean, test_std], test_set.targets, batch_size, verbose=0))
Out[24]:
We get very similar results to the theano implementation.
In [28]:
model_keras = build_EMOS_network_keras(compile=True, optimizer='sgd', lr=0.1)
In [29]:
date_str_start = '2016-01-01'
date_str_stop = '2017-01-01'
# This function loops over the days.
train_crps_list, valid_crps_list, results_df = loop_over_days(
DATA_DIR,
model_keras,
date_str_start, date_str_stop,
window_size=window_size,
fclt=fclt,
epochs_max=steps_max,
early_stopping_delta=early_stopping_delta,
lr=lr,
verbose=0,
model_type='EMOS_Network_keras')
The keras implementation is slower than the pure theano version. This could be due to the overhead of calling model.fit many many times.
In [30]:
np.mean(train_crps_list), np.mean(valid_crps_list)
Out[30]:
The results are slightly better than the theano implementation. This could be due to random variability or due to the difference in the early stopping algorithm.
In [31]:
results_df.to_csv(results_dir + 'emos_network_rolling_window_keras.csv')
In [32]:
train_dates = ['2015-01-01', '2016-01-01']
test_dates = ['2016-01-01', '2017-01-01']
In [33]:
# Load data sets
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates, test_dates)
In [34]:
model_keras = build_EMOS_network_keras(compile=True, optimizer='adam', lr=0.1)
In [35]:
# Split dataset
train_mean = train_set.features[:, 0]
train_std = train_set.features[:, 1]
test_mean = test_set.features[:, 0]
test_std = test_set.features[:, 1]
In [36]:
model_keras.fit([train_mean, train_std], train_set.targets, epochs=10,
batch_size=1024,
validation_data=[[test_mean, test_std], test_set.targets])
Out[36]:
So we get a very similar CRPS compared to the 25 day rolling window. This suggests that the seasonality is not that important.
In [37]:
# Get predictions
preds = model_keras.predict([test_mean, test_std])
In [38]:
# Save predictions
results_df = create_results_df(test_set.date_strs, test_set.station_ids,
preds[:, 0], preds[:, 1])
In [39]:
results_df.head()
Out[39]:
In [40]:
results_df.to_csv(results_dir + 'emos_network_train_2015_pred_2016.csv')
In [ ]: