<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 # Table of Contents
<div class="toc" style="margin-top
Step1:  Set up data
This is now also done inside the get_train_test_sets function. seq_len is the number of timesteps (including the one to predict). We will start out with a moderate length of 5 days, training for 2015, predicting for 2016.
Step2:  The arrays have dimensions [sample, time step, feature]
RNNs with only temperature data
As a comparison. Our simple networks got a train/test loss of around 1.07/1.01.
I am using a Gated Recurrent Unit (GRU) as my recurrent layer. LSTM is probably the more common one, but GRU is slightly cheaper and for our simple applications provides similar results. 
RNN predicting only the last target
Step4:  So we get a better train score and a worse validation score. This indicates overfitting. 
Sequence RNN
Step5:  Same as with the first RNN above we seem to overfit to the dataset, but maybe not as strongly. Let's now try a more complex model with a longer sequence length.
Longer sequence
Step6:  So again we are overfitting, but maybe there is something to be learned. Let's first add some regularization and then try a longer training set.
Add dropout
Step7:  So with drop out we get slightly better validation results, but we are still starting to overfit. I think there is a lot of parameter tuning that would be possible with the complexity of the network and so forth. 
Reference experiment with longer training set.
Step8:  Maybe a small improvement. Now let's test our sequence model with a longer training period.
Sequence model with longer training set
Step9:  Get additional variables
Step10:  Combining with embeddings
Step11:  Embedding with only temperature
Step12:  Add previous forecasts and observations as features
Let's try a very simple form of memory. Simply adding the forecast, observation and error at the current time.
Step13:  Erm ok wow, that is pretty incredible. But wait maybe this is very similar to the embedding information.
Step14:  Alright, I need to check whether I am cheating, but for now let's try to build the best model.
Step15:  Go further back and test importance
<ASSISTANT_TASK:>
Python Code:
# Imports
import sys
sys.path.append('../')   # This is where all the python files are!
from importlib import reload
import utils; reload(utils)
from utils import *
import keras_models; reload(keras_models)
from keras_models import *
import losses; reload(losses)
from losses import crps_cost_function, crps_cost_function_seq
import matplotlib.pyplot as plt
%matplotlib inline
import keras
from keras.layers import Input, Dense, merge, Embedding, Flatten, Dropout, \
    SimpleRNN, LSTM, TimeDistributed, GRU, Dropout, Masking
from keras.layers.merge import Concatenate
from keras.models import Model, Sequential
import keras.backend as K
from keras.callbacks import EarlyStopping
from keras.optimizers import SGD, Adam
# Use this if you want to limit the GPU RAM usage
import tensorflow as tf
from keras.backend.tensorflow_backend import set_session
config = tf.ConfigProto()
config.gpu_options.allow_growth = True
set_session(tf.Session(config=config))
# Basic setup
DATA_DIR = '/Volumes/STICK/data/ppnn_data/'  # Mac
# DATA_DIR = '/project/meteo/w2w/C7/ppnn_data/'   # LMU
results_dir = '../results/'
window_size = 25   # Days in rolling window
fclt = 48   # Forecast lead time in hours
Explanation: # Table of Contents
<div class="toc" style="margin-top: 1em;"><ul class="toc-item" id="toc-level0"><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Set-up-data" data-toc-modified-id="Set-up-data-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Set up data</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#RNNs-with-only-temperature-data" data-toc-modified-id="RNNs-with-only-temperature-data-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>RNNs with only temperature data</a></span><ul class="toc-item"><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#RNN-predicting-only-the-last-target" data-toc-modified-id="RNN-predicting-only-the-last-target-2.1"><span class="toc-item-num">2.1&nbsp;&nbsp;</span>RNN predicting only the last target</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Sequence-RNN" data-toc-modified-id="Sequence-RNN-2.2"><span class="toc-item-num">2.2&nbsp;&nbsp;</span>Sequence RNN</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Longer-sequence" data-toc-modified-id="Longer-sequence-2.3"><span class="toc-item-num">2.3&nbsp;&nbsp;</span>Longer sequence</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Add-dropout" data-toc-modified-id="Add-dropout-2.4"><span class="toc-item-num">2.4&nbsp;&nbsp;</span>Add dropout</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Reference-experiment-with-longer-training-set." data-toc-modified-id="Reference-experiment-with-longer-training-set.-2.5"><span class="toc-item-num">2.5&nbsp;&nbsp;</span>Reference experiment with longer training set.</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Sequence-model-with-longer-training-set" data-toc-modified-id="Sequence-model-with-longer-training-set-2.6"><span class="toc-item-num">2.6&nbsp;&nbsp;</span>Sequence model with longer training set</a></span></li></ul></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Get-additional-variables" data-toc-modified-id="Get-additional-variables-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>Get additional variables</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Combining-with-embeddings" data-toc-modified-id="Combining-with-embeddings-4"><span class="toc-item-num">4&nbsp;&nbsp;</span>Combining with embeddings</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Embedding-with-only-temperature" data-toc-modified-id="Embedding-with-only-temperature-5"><span class="toc-item-num">5&nbsp;&nbsp;</span>Embedding with only temperature</a></span></li><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Add-previous-forecasts-and-observations-as-features" data-toc-modified-id="Add-previous-forecasts-and-observations-as-features-6"><span class="toc-item-num">6&nbsp;&nbsp;</span>Add previous forecasts and observations as features</a></span><ul class="toc-item"><li><span><a href="http://localhost:8888/notebooks/notebooks/rnn.ipynb#Go-further-back-and-test-importance" data-toc-modified-id="Go-further-back-and-test-importance-6.1"><span class="toc-item-num">6.1&nbsp;&nbsp;</span>Go further back and test importance</a></span></li></ul></li></ul></div>
Recurrent neural networks
In this notebook we will try out RNNs for our post-processing. The idea here is that there might be some extra information in looking at data from previous time steps.
RNNs take quite a long time to train, so I am using a GPU here.
End of explanation
seq_len=5
train_dates = ['2015-01-01', '2016-01-01']
test_dates =  ['2016-01-01', '2017-01-01']
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates, test_dates, 
                                          seq_len=seq_len, fill_value=-999.)
train_set.features.shape, train_set.targets.shape
Explanation: Set up data
This is now also done inside the get_train_test_sets function. seq_len is the number of timesteps (including the one to predict). We will start out with a moderate length of 5 days, training for 2015, predicting for 2016.
End of explanation
batch_size = 1024
hidden_nodes = 100   # Number of hidden nodes inside RNN cell
inp = Input(shape=(seq_len, 2, )) # time step, feature
x = GRU(hidden_nodes)(inp)
x = Dense(2, activation='linear')(x)
rnn_model = Model(inputs=inp, outputs=x)
rnn_model.compile(optimizer=Adam(0.01), loss=crps_cost_function)
rnn_model.summary()
rnn_model.fit(train_set.features, train_set.targets[:,-1], epochs=10, batch_size=batch_size,
              validation_data=(test_set.features, test_set.targets[:,-1]))
Explanation: The arrays have dimensions [sample, time step, feature]
RNNs with only temperature data
As a comparison. Our simple networks got a train/test loss of around 1.07/1.01.
I am using a Gated Recurrent Unit (GRU) as my recurrent layer. LSTM is probably the more common one, but GRU is slightly cheaper and for our simple applications provides similar results. 
RNN predicting only the last target
End of explanation
inp = Input(shape=(seq_len, 2, )) # time step, feature
x = GRU(hidden_nodes, return_sequences=True)(inp)
x = TimeDistributed(Dense(2, activation='linear'))(x)
seq_rnn_model = Model(inputs=inp, outputs=x)
seq_rnn_model.summary()
seq_rnn_model.compile(optimizer=Adam(0.01), loss=crps_cost_function_seq, 
                      sample_weight_mode="temporal")
def train_and_valid(model, train_set, test_set, epochs, batch_size, verbose=0, emb=False):
    Write our own function to train and validate, 
    because the keras fit function cannot handle sample weights for training
    and validation at the same time.
    
    train_inp = [train_set.features, train_set.cont_ids] if emb else train_set.features
    test_inp = [test_set.features, test_set.cont_ids] if emb else test_set.features
    for i in range(epochs):
        print('Epoch:', i+1)
        t1 = timeit.default_timer()
        
        h = model.fit(train_inp, train_set.targets, epochs=1, batch_size=batch_size, 
                      sample_weight=train_set.sample_weights, verbose=verbose)
        t2 = timeit.default_timer()
        print('Train loss: %.4f - Valid loss: %.4f - Time: %.1fs' % (h.history['loss'][0],  
                    model.evaluate(test_inp, test_set.targets, batch_size=10000, 
                       sample_weight=test_set.sample_weights, verbose=verbose), 
                    t2 - t1))
train_and_valid(seq_rnn_model, train_set, test_set, 10, batch_size)
Explanation: So we get a better train score and a worse validation score. This indicates overfitting. 
Sequence RNN
End of explanation
seq_len = 20
train_dates = ['2015-01-01', '2016-01-01']
test_dates =  ['2016-01-01', '2017-01-01']
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates, test_dates, 
                                          seq_len=seq_len, fill_value=-999.)
hidden_nodes = 200
inp = Input(shape=(seq_len, 2, )) # time step, feature
x = GRU(hidden_nodes, return_sequences=True)(inp)
x = TimeDistributed(Dense(2, activation='linear'))(x)
seq_rnn_model = Model(inputs=inp, outputs=x)
seq_rnn_model.summary()
seq_rnn_model.compile(optimizer=Adam(0.01), loss=crps_cost_function_seq, 
                      sample_weight_mode="temporal")
train_and_valid(seq_rnn_model, train_set, test_set, 10, batch_size)
Explanation: Same as with the first RNN above we seem to overfit to the dataset, but maybe not as strongly. Let's now try a more complex model with a longer sequence length.
Longer sequence
End of explanation
inp = Input(shape=(seq_len, 2, )) # time step, feature
x = GRU(hidden_nodes, return_sequences=True, recurrent_dropout=0.5)(inp)
x = TimeDistributed(Dense(2, activation='linear'))(x)
seq_rnn_model = Model(inputs=inp, outputs=x)
seq_rnn_model.compile(optimizer=Adam(0.001), loss=crps_cost_function_seq, 
                      sample_weight_mode="temporal")
train_and_valid(seq_rnn_model, train_set, test_set, 10, batch_size)
Explanation: So again we are overfitting, but maybe there is something to be learned. Let's first add some regularization and then try a longer training set.
Add dropout
End of explanation
train_dates_long = ['2008-01-01', '2016-01-01']
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates_long, test_dates)
# Copied from fc_network notebook
def build_fc_model():
    inp = Input(shape=(2,))
    x = Dense(2, activation='linear')(inp)
    return Model(inputs=inp, outputs=x)
fc_model = build_fc_model()
fc_model.compile(optimizer=Adam(0.1), loss=crps_cost_function)
fc_model.fit(train_set.features, train_set.targets, epochs=10, batch_size=1024,
             validation_data=[test_set.features, test_set.targets])
Explanation: So with drop out we get slightly better validation results, but we are still starting to overfit. I think there is a lot of parameter tuning that would be possible with the complexity of the network and so forth. 
Reference experiment with longer training set.
End of explanation
seq_len = 20
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates_long, test_dates, 
                                          seq_len=seq_len, fill_value=-999.)
def build_seq_rnn(hidden_nodes, n_features, dropout=0, lr=0.01):
    inp = Input(shape=(seq_len, n_features, )) # time step, feature
    x = GRU(hidden_nodes, return_sequences=True, recurrent_dropout=dropout)(inp)
    x = TimeDistributed(Dense(2, activation='linear'))(x)
    seq_rnn_model = Model(inputs=inp, outputs=x)
    seq_rnn_model.compile(optimizer=Adam(lr), loss=crps_cost_function_seq, 
                          sample_weight_mode="temporal")
    return seq_rnn_model
inp = Input(shape=(seq_len, 2, )) # time step, feature
x = GRU(hidden_nodes, return_sequences=True, recurrent_dropout=0.5)(inp)
x = TimeDistributed(Dense(2, activation='linear'))(x)
seq_rnn_model = Model(inputs=inp, outputs=x)
seq_rnn_model.compile(optimizer=Adam(0.001), loss=crps_cost_function_seq, 
                      sample_weight_mode="temporal")
# This takes several minutes on the GPU
# Epoch counter: 7
train_and_valid(seq_rnn_model, train_set, test_set, 2, batch_size)
Explanation: Maybe a small improvement. Now let's test our sequence model with a longer training period.
Sequence model with longer training set
End of explanation
from collections import OrderedDict
aux_dict = OrderedDict()
aux_dict['data_aux_geo_interpolated.nc'] = ['orog', 
                                            'station_alt', 
                                            'station_lat', 
                                            'station_lon']
aux_dict['data_aux_pl500_interpolated_00UTC.nc'] = ['u_pl500_fc',
                                                    'v_pl500_fc',
                                                    'gh_pl500_fc']
aux_dict['data_aux_pl850_interpolated_00UTC.nc'] = ['u_pl850_fc',
                                                    'v_pl850_fc',
                                                    'q_pl850_fc']
aux_dict['data_aux_surface_interpolated_00UTC.nc'] = ['cape_fc',
                                                      'sp_fc',
                                                      'tcc_fc']
# Start with just one training year
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates, test_dates, 
                                          seq_len=seq_len, fill_value=-999., aux_dict=aux_dict)
train_set.cont_ids.shape
n_features = train_set.features.shape[-1]
n_features
seq_rnn_model = build_seq_rnn(hidden_nodes, n_features, dropout=0.5, lr=0.01)
# Epoch counter: 8
train_and_valid(seq_rnn_model, train_set, test_set, 2, batch_size, verbose=1)
seq_rnn_model.summary()
Explanation: Get additional variables
End of explanation
def build_seq_rnn_with_embeddings(seq_len, hidden_nodes, n_features, emb_size, max_id, 
                                  recurrent_dropout=0, dropout=0, lr=0.01):
    features_inp = Input(shape=(seq_len, n_features, )) # time step, feature
    id_in = Input(shape=(seq_len,))
    emb = Embedding(max_id + 1, emb_size)(id_in)
    x = GRU(hidden_nodes, return_sequences=True, recurrent_dropout=recurrent_dropout)(features_inp)
    x = Concatenate()([x, emb])
    x = Dropout(dropout)(x)
    x = TimeDistributed(Dense(2, activation='linear'))(x)
    model = Model(inputs=[features_inp, id_in], outputs=x)
    model.compile(optimizer=Adam(lr), loss=crps_cost_function_seq, 
                          sample_weight_mode="temporal")
    return model
emb_size = 5
max_id = int(np.max([train_set.cont_ids.max(), test_set.cont_ids.max()]))
hidden_nodes, max_id
emb_rnn = build_seq_rnn_with_embeddings(hidden_nodes, n_features, emb_size, max_id, 0.5)
emb_rnn.summary()
train_and_valid(emb_rnn, train_set, test_set, 5, batch_size, emb=True)
Explanation: Combining with embeddings
End of explanation
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates, test_dates, 
                                          seq_len=2, fill_value=-999.)
emb_rnn = build_seq_rnn_with_embeddings(2, 10, 2, emb_size, max_id, recurrent_dropout=0)
emb_rnn.summary()
train_and_valid(emb_rnn, train_set, test_set, 5, batch_size, emb=True)
train_dates_long = ['2008-01-01', '2016-01-01']
train_set_long, test_set = get_train_test_sets(DATA_DIR, train_dates_long, test_dates, 
                                          seq_len=seq_len, fill_value=-999.)
emb_rnn = build_seq_rnn_with_embeddings(seq_len, 30, 2, emb_size, max_id, recurrent_dropout=0.3)
train_and_valid(emb_rnn, train_set_long, test_set, 15, batch_size, emb=True)
Explanation: Embedding with only temperature
End of explanation
train_dates = ['2015-01-01', '2016-01-01']
test_dates =  ['2016-01-01', '2017-01-01']
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates, test_dates,
                                          add_current_error=True)
train_set.feature_names
train_set.features.shape
test_set.features[10000:11000, :]
fc_model = build_fc_model(5, 2, compile=True)
# Note: I am running this cell several times 
fc_model.fit(train_set.features, train_set.targets, epochs=10, batch_size=1024,
             validation_data=[test_set.features, test_set.targets])
Explanation: Add previous forecasts and observations as features
Let's try a very simple form of memory. Simply adding the forecast, observation and error at the current time.
End of explanation
emb_size = 3
max_id = int(np.max([train_set.cont_ids.max(), test_set.cont_ids.max()]))
max_id
emb_model = build_emb_model(5, 2, [], emb_size, max_id, compile=True,
                            lr=0.01)
emb_model.fit([train_set.features, train_set.cont_ids], train_set.targets, 
              epochs=10, batch_size=1024, 
              validation_data=[[test_set.features, test_set.cont_ids], test_set.targets])
Explanation: Erm ok wow, that is pretty incredible. But wait maybe this is very similar to the embedding information.
End of explanation
from collections import OrderedDict
more_aux_dict = OrderedDict()
more_aux_dict['data_aux_geo_interpolated.nc'] = ['orog', 
                                            'station_alt', 
                                            'station_lat', 
                                            'station_lon']
more_aux_dict['data_aux_pl500_interpolated_00UTC.nc'] = ['u_pl500_fc',
                                                    'v_pl500_fc',
                                                    'gh_pl500_fc']
more_aux_dict['data_aux_pl850_interpolated_00UTC.nc'] = ['u_pl850_fc',
                                                    'v_pl850_fc',
                                                    'q_pl850_fc']
more_aux_dict['data_aux_surface_interpolated_00UTC.nc'] = ['cape_fc',
                                                      'sp_fc',
                                                      'tcc_fc']
more_aux_dict['data_aux_surface_more_interpolated_part1_00UTC.nc']  = [
    'sshf_fc', 'slhf_fc', 'u10_fc','v10_fc'
]
more_aux_dict['data_aux_surface_more_interpolated_part2_00UTC.nc']  = [
    'ssr_fc', 'str_fc', 'd2m_fc','sm_fc'
]
train_dates = ['2015-01-01', '2016-01-01']
test_dates =  ['2016-01-01', '2017-01-01']
more_train_set, more_test_set = get_train_test_sets(DATA_DIR, train_dates, test_dates,
                                         aux_dict=more_aux_dict, add_current_error=True)
emb_size = 3
max_id = int(np.max([more_train_set.cont_ids.max(), more_test_set.cont_ids.max()]))
max_id
emb_model = build_emb_model(more_train_set.features.shape[1], 2, [50], 3, max_id, 
                            compile=True, lr=0.01)
from IPython.display import SVG, 
from keras.utils.vis_utils import model_to_dot
SVG(model_to_dot(emb_model, show_shapes=True).create(prog='dot', format='svg'))
emb_model.fit([more_train_set.features, more_train_set.cont_ids], more_train_set.targets, epochs=30, 
              batch_size=4096, validation_split=0.0)
              #callbacks=[EarlyStopping(monitor='val_loss', 
              #                         min_delta=0,
              #                         patience=3)])
emb_model.evaluate([more_test_set.features, more_test_set.cont_ids], more_test_set.targets, batch_size=10000)
emb_model.summary()
long_train_dates = ['2008-01-01', '2016-01-01']
long_more_train_set, more_test_set = get_train_test_sets(DATA_DIR, long_train_dates, test_dates,
                                                         aux_dict=more_aux_dict,
                                                        add_current_error=True)
emb_model = build_emb_model(long_more_train_set.features.shape[1], 2, [50], 3, max_id, 
                            compile=True, lr=0.01)
emb_model.fit([long_more_train_set.features, long_more_train_set.cont_ids], long_more_train_set.targets, epochs=50, 
              batch_size=4096, validation_split=0.2,
              callbacks=[EarlyStopping(monitor='val_loss', 
                                       min_delta=0,
                                       patience=2)])
emb_model.evaluate([more_test_set.features, more_test_set.cont_ids], test_set.targets, batch_size=10000)
Explanation: Alright, I need to check whether I am cheating, but for now let's try to build the best model.
End of explanation
# Test current error
train_dates = ['2015-01-01', '2016-01-01']
test_dates =  ['2016-01-01', '2017-01-01']
train_set, test_set = get_train_test_sets(DATA_DIR, train_dates, test_dates,
                                          add_current_error=True, 
                                          current_error_len=1)
train_set.feature_names
fc_model = build_fc_model(4, 2, compile=True)
fc_model.fit(train_set.features, train_set.targets, epochs=10, batch_size=1024,
             validation_data=[test_set.features, test_set.targets])
Explanation: Go further back and test importance
End of explanation
<END_TASK>