import numpy as np
import keras
import tensorflow as tf
from matplotlib import pyplot as plt
import pandas as pd
from keras.layers.embeddings import Embedding
from keras.layers import concatenate, Lambda
import os, sys
from weather_model import Seq2Seq_MVE_subnets_swish, weather_conv1D, CausalCNN, RNN_builder, Seq2Seq, Seq2Seq_MVE, Seq2Seq_MVE_subnets
from keras.models import load_model, model_from_json
#from utils import random_sine, plot_prediction
#learning_rate = 0.01
#decay = 0 # Learning rate decay
model_save_path = '../models/'
#loss = "mse" # Other loss functions are possible, see Keras documentation.
# Regularisation isn't really needed for this application
#lambda_regulariser = None #0.000001 # Will not be used if regulariser is None
#regulariser = None # Possible regulariser: keras.regularizers.l2(lambda_regulariser)
#steps_per_epoch = 200 # batch_size * steps_per_epoch = total number of training examples
#num_signals = 2 # The number of random sine waves the compose the signal. The more sine waves, the harder the problem.
def crop(dimension, start, end):
# Crops (or slices) a Tensor on a given dimension from start to end
# example : to crop tensor x[:, :, 5:10]
# call slice(2, 5, 10) as you want to crop on the second dimension
def func(x):
if dimension == 0:
return x[start: end]
if dimension == 1:
return x[:, start: end]
if dimension == 2:
return x[:, :, start: end]
if dimension == 3:
return x[:, :, :, start: end]
if dimension == 4:
return x[:, :, :, :, start: end]
return Lambda(func)
class WeatherConv1D:
def __init__(self, regulariser=None,
lr=0.001, decay=0, loss="mse",
layers=[35, 35], batch_size=256,
input_len=37, input_features=29,
strides_len=3, kernel_size=5):
self.regulariser = regulariser
self.layers = layers
self.lr = lr
self.decay = decay
self.loss = loss
self.pred_result = None
self.batch_size = batch_size
self.input_len = input_len
self.input_features = input_features
self.kernel_strides = strides_len
self.kernel_size = kernel_size
print('Initialized!')
def build_graph(self):
keras.backend.clear_session() # clear session/graph
self.model = weather_conv1D(self.layers, self.lr,
self.decay, self.loss, self.input_len,
self.input_features, self.kernel_strides, self.kernel_size)
print(self.model.summary())
def sample_batch(self, data_inputs, ground_truth, ruitu_inputs, batch_size, certain_id=None, certain_feature=None):
max_i, _, max_j, _ = data_inputs.shape # Example: (1148, 37, 10, 9)-(sample_ind, timestep, sta_id, features)
if certain_id == None and certain_feature == None:
id_ = np.random.randint(max_j, size=batch_size)
i = np.random.randint(max_i, size=batch_size)
batch_ouputs = ground_truth[i,:,id_,:]
batch_ruitu = ruitu_inputs[i,:,id_,:]
elif certain_id != None:
pass
return batch_ruitu, batch_ouputs
def order_batch(self, data_inputs, ground_truth, ruitu_inputs, batch_size, certain_id=None, certain_feature=None):
pass #TODO:
def fit(self, train_input_ruitu, train_labels,
val_input_ruitu, val_labels, batch_size,
iterations=300, validation=True):
self.optimizer = keras.optimizers.Adam(lr=self.lr, decay=self.decay)
self.model.compile(optimizer = self.optimizer, loss=self.loss)
print('Train batch size: {}'.format(batch_size))
print('Validation on data size of {};'.format(val_input_ruitu.shape[0]))
for i in range(iterations):
batch_ruitu, batch_labels = self.sample_batch(train_input_ruitu, train_labels,
train_input_ruitu, batch_size=batch_size)
loss_ = self.model.train_on_batch(x=[batch_ruitu],
y=[batch_labels])
if (i+1)%50 == 0:
print('Iteration:{}/{}. Training batch loss:{}'.
format(i+1, iterations, loss_))
if validation :
self.evaluate(val_input_ruitu, val_labels, each_station_display=False)
print('###'*10)
print('Train finish! Total validation loss:')
self.evaluate(val_input_ruitu, val_labels, each_station_display=True)
def evaluate(self, data_input_ruitu, data_labels, each_station_display=False):
all_loss=[]
for i in range(10): # iterate for each station. (sample_ind, timestep, staionID, features)
val_loss= self.model.evaluate(x=[data_input_ruitu[:,:,i,:]],
y=[data_labels[:,:,i,:]], verbose=False)
all_loss.append(val_loss)
if each_station_display:
print('\tFor station 9000{}, evaluated loss: {}'.format(i+1, val_loss))
print('Mean evaluated loss on all stations:', np.mean(all_loss))
#return np.mean(all_loss)
def predict(self, batch_ruitu):
pred_result_list = []
for i in range(10):
#print('Predict for station: 9000{}'.format(i+1))
result = self.model.predict(x=[batch_ruitu[:,:,i,:]])
result = np.squeeze(result, axis=0)
#all_pred[i] = result
pred_result_list.append(result)
#pass
pred_result = np.stack(pred_result_list, axis=0)
#return all_pred, pred_result
print('Predict shape (10,37,3) means (stationID, timestep, features). Features include: t2m, rh2m and w10m')
self.pred_result = pred_result
return pred_result
def renorm_for_submit(self, pred_mean, pred_var):
if self.pred_result is None:
print('You must run self.predict(batch_inputs, batch_ruitu) firstly!!')
else:
df_empty = pd.DataFrame(columns=['FORE_data', 't2m', 'rh2m', 'w10m'])
target_list=['t2m','rh2m','w10m']
self.obs_range_dic={'t2m':[-30,42], # Official value: [-20,42]
'rh2m':[0.0,100.0],
'w10m':[0.0, 30.0]}
for j, target_v in enumerate(self.target_list):
series_ids = pd.Series()
series_targets = pd.Series()
renorm_value = renorm(self.pred_result[:,:,j], self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
for i in range(10):
if i != 9:
id_num = '0'+str(i+1)
else:
id_num = str(10)
sta_name_time = '900'+id_num+'_'
time_str_list=[]
for t in range(37):
if t < 10:
time_str= sta_name_time + '0'+ str(t)
else:
time_str = sta_name_time + str(t)
time_str_list.append(time_str)
series_id = pd.Series(time_str_list)
series_target = pd.Series(renorm_value[i])
series_ids = pd.concat([series_ids, series_id])
series_targets = pd.concat([series_targets, series_target])
df_empty['FORE_data'] = series_ids
df_empty[target_v] = series_targets
return df_empty
class CausalCNN_Class(WeatherConv1D):
def __init__(self, regulariser,lr, decay, loss,
n_filters, strides_len, kernel_size, seq_len,
input_features, output_features, dilation_rates):
self.regulariser=regulariser
self.n_filters=n_filters
self.lr=lr
self.decay=decay
self.loss=loss
self.seq_len=seq_len
self.input_features=input_features
self.output_features = output_features
self.strides_len=strides_len
self.kernel_size=kernel_size
self.dilation_rates=dilation_rates
def build_graph(self):
keras.backend.clear_session() # clear session/graph
self.model = CausalCNN(self.n_filters, self.lr,
self.decay, self.loss,
self.seq_len, self.input_features,
self.strides_len, self.kernel_size,
self.dilation_rates)
print(self.model.summary())
class FNN(WeatherConv1D):
def __init__(self, regulariser,lr, decay, loss,
layers, batch_size, seq_len, input_features, output_features):
self.regulariser=regulariser
self.layers=layers
self.lr=lr
self.decay=decay
self.loss=loss
self.seq_len=seq_len
self.input_features=input_features
self.output_features = output_features
def build_graph(self):
keras.backend.clear_session() # clear session/graph
self.model = weather_fnn(self.layers, self.lr,
self.decay, self.loss, self.seq_len,
self.input_features, self.output_features)
print(self.model.summary())
class Enc_Dec:
def __init__(self, num_input_features, num_output_features, num_decoder_features,
input_sequence_length, target_sequence_length,
num_steps_to_predict, regulariser = None,
lr=0.001, decay=0, loss = "mse",
layers=[35, 35]):
self.num_input_features = num_input_features
self.num_output_features = num_output_features
self.num_decoder_features = num_decoder_features
self.input_sequence_length = input_sequence_length
self.target_sequence_length = target_sequence_length
self.num_steps_to_predict = num_steps_to_predict
self.regulariser = regulariser
self.layers = layers
self.lr = lr
self.decay = decay
self.loss = loss
self.pred_result = None
self.train_loss=[]
self.target_list=['t2m','rh2m','w10m']
self.obs_range_dic={'t2m':[-30,42], # Official value: [-20,42]
'rh2m':[0.0,100.0],
'w10m':[0.0, 30.0]}
print('Initialized!')
def build_graph(self):
keras.backend.clear_session() # clear session/graph
self.optimiser = keras.optimizers.Adam(lr=self.lr, decay=self.decay)
# Define an input sequence.
encoder_inputs = keras.layers.Input(shape=(None, self.num_input_features), name='encoder_inputs')
# Create a list of RNN Cells, these are then concatenated into a single layer
# with the RNN layer.
encoder_cells = []
for hidden_neurons in self.layers:
encoder_cells.append(keras.layers.GRUCell(hidden_neurons,
kernel_regularizer = self.regulariser,
recurrent_regularizer = self.regulariser,
bias_regularizer = self.regulariser))
encoder = keras.layers.RNN(encoder_cells, return_state=True)
encoder_outputs_and_states = encoder(encoder_inputs)
# Discard encoder outputs and only keep the states.
encoder_states = encoder_outputs_and_states[1:]
# Define a decoder sequence.
decoder_inputs = keras.layers.Input(shape=(None, self.num_decoder_features), name='decoder_inputs')
decoder_cells = []
for hidden_neurons in self.layers:
decoder_cells.append(keras.layers.GRUCell(hidden_neurons,
kernel_regularizer = self.regulariser,
recurrent_regularizer = self.regulariser,
bias_regularizer = self.regulariser))
decoder = keras.layers.RNN(decoder_cells, return_sequences=True, return_state=True)
# Set the initial state of the decoder to be the ouput state of the encoder.
decoder_outputs_and_states = decoder(decoder_inputs, initial_state=encoder_states)
# Only select the output of the decoder (not the states)
decoder_outputs = decoder_outputs_and_states[0]
# Apply a dense layer with linear activation to set output to correct dimension
# and scale (tanh is default activation for GRU in Keras, our output sine function can be larger then 1)
decoder_dense1 = keras.layers.Dense(units=64,
activation='tanh',
kernel_regularizer = self.regulariser,
bias_regularizer = self.regulariser, name='dense_tanh')
output_dense = keras.layers.Dense(self.num_output_features,
activation='sigmoid',
kernel_regularizer = self.regulariser,
bias_regularizer = self.regulariser, name='output_sig')
#densen1=decoder_dense1(decoder_outputs)
decoder_outputs = output_dense(decoder_outputs)
# Create a model using the functional API provided by Keras.
self.model = keras.models.Model(inputs=[encoder_inputs, decoder_inputs], outputs=decoder_outputs)
print(self.model.summary())
def sample_batch(self, data_inputs, ground_truth, ruitu_inputs, batch_size, certain_id=None, certain_feature=None):
max_i, _, max_j, _ = data_inputs.shape # Example: (1148, 37, 10, 9)-(sample_ind, timestep, sta_id, features)
if certain_id == None and certain_feature == None:
id_ = np.random.randint(max_j, size=batch_size)
i = np.random.randint(max_i, size=batch_size)
batch_inputs = data_inputs[i,:,id_,:]
batch_ouputs = ground_truth[i,:,id_,:]
batch_ruitu = ruitu_inputs[i,:,id_,:]
elif certain_id != None:
pass
return batch_inputs, batch_ruitu, batch_ouputs
def fit(self, train_input_obs, train_input_ruitu, train_labels,
val_input_obs, val_input_ruitu, val_labels, batch_size,
iterations=300, validation=True):
self.model.compile(optimizer = self.optimiser, loss=self.loss)
print('Train batch size: {}'.format(batch_size))
print('Validation on data size of {};'.format(val_input_obs.shape[0]))
for i in range(iterations):
batch_inputs, batch_ruitu, batch_labels = self.sample_batch(train_input_obs, train_labels,
train_input_ruitu, batch_size=batch_size)
loss_ = self.model.train_on_batch(x=[batch_inputs, batch_ruitu],
y=[batch_labels])
if (i+1)%50 == 0:
print('Iteration:{}/{}. Training batch loss:{}'.
format(i+1, iterations, loss_))
if validation :
self.evaluate(val_input_obs, val_input_ruitu, val_labels, each_station_display=False)
print('###'*10)
print('Train finish! Total validation loss:')
self.evaluate(val_input_obs, val_input_ruitu, val_labels, each_station_display=True)
def evaluate(self, data_input_obs, data_input_ruitu, data_labels, each_station_display=False):
assert data_input_ruitu.shape[0] == data_input_obs.shape[0] == data_labels.shape[0], 'Shape Error'
#assert data_input_obs.shape[1] == 28 and data_input_obs.shape[2] == 10 and data_input_obs.shape[3] == 9, 'Error! Obs input shape must be (None, 28,10,9)'
assert data_input_ruitu.shape[1] == 37 and data_input_ruitu.shape[2] == 10 and data_input_ruitu.shape[3] == 29, 'Error! Ruitu input shape must be (None, 37,10,29)'
assert data_labels.shape[1] == 37 and data_labels.shape[2] == 10 and data_labels.shape[3] == 3, 'Error! Ruitu input shape must be (None, 37,10,3)'
all_loss=[]
for i in range(10): # iterate for each station. (sample_ind, timestep, staionID, features)
val_loss= self.model.evaluate(x=[data_input_obs[:,:,i,:], data_input_ruitu[:,:,i,:]],
y=[data_labels[:,:,i,:]], verbose=False)
all_loss.append(val_loss)
if each_station_display:
print('\tFor station 9000{}, evaluated loss: {}'.format(i+1, val_loss))
print('Mean evaluated loss on all stations:', np.mean(all_loss))
def predict(self, batch_inputs, batch_ruitu):
assert batch_ruitu.shape[0] == batch_inputs.shape[0], 'Shape Error'
assert batch_inputs.shape[1] == 28 and batch_inputs.shape[2] == 10 and batch_inputs.shape[3] == 9, 'Error! Obs input shape must be (None, 28,10,9)'
assert batch_ruitu.shape[1] == 37 and batch_ruitu.shape[2] == 10 and batch_ruitu.shape[3] == 29, 'Error! Ruitu input shape must be (None, 37,10, 29)'
#all_pred={}
pred_result_list = []
for i in range(10):
#print('Predict for station: 9000{}'.format(i+1))
result = self.model.predict(x=[batch_inputs[:,:,i,:], batch_ruitu[:,:,i,:]])
result = np.squeeze(result, axis=0)
#all_pred[i] = result
pred_result_list.append(result)
#pass
pred_result = np.stack(pred_result_list, axis=0)
#return all_pred, pred_result
print('Predict shape (10,37,3) means (stationID, timestep, features). Features include: t2m, rh2m and w10m')
self.pred_result = pred_result
return pred_result
def renorm_for_submit(self, pred_mean, pred_var=None):
'''
# TODO: Add three strategies for output
'''
assert self.pred_result is not None, 'You must run self.predict(batch_inputs, batch_ruitu) firstly!!'
assert pred_mean.shape == (10, 37, 3), 'Error! This funtion ONLY works for one data sample with shape (10, 37, 3). Any data shape (None, 10, 37, 3) will leads this error!'
df_empty = pd.DataFrame(columns=['FORE_data', 't2m', 'rh2m', 'w10m'])
for j, target_v in enumerate(self.target_list):
series_ids = pd.Series()
series_targets = pd.Series()
renorm_value = renorm(pred_mean[:,:,j], self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
for i in range(10):
if i != 9:
id_num = '0'+str(i+1)
else:
id_num = str(10)
sta_name_time = '900'+id_num+'_'
time_str_list=[]
for t in range(37):
if t < 10:
time_str= sta_name_time + '0'+ str(t)
else:
time_str = sta_name_time + str(t)
time_str_list.append(time_str)
series_id = pd.Series(time_str_list)
series_target = pd.Series(renorm_value[i])
series_ids = pd.concat([series_ids, series_id])
series_targets = pd.concat([series_targets, series_target])
df_empty['FORE_data'] = series_ids
df_empty[target_v] = series_targets
return df_empty
#pass
def plot_prediction(self, x, y_true, y_pred, input_ruitu=None):
"""Plots the predictions.
Arguments
---------
x: Input sequence of shape (input_sequence_length,
dimension_of_signal)
y_true: True output sequence of shape (input_sequence_length,
dimension_of_signal)
y_pred: Predicted output sequence (input_sequence_length,
dimension_of_signal)
input_ruitu: Ruitu output sequence
"""
plt.figure(figsize=(12, 3))
output_dim = x.shape[-1]# feature dimension
for j in range(output_dim):
past = x[:, j]
true = y_true[:, j]
pred = y_pred[:, j]
if input_ruitu is not None:
ruitu = input_ruitu[:, j]
label1 = "Seen (past) values" if j==0 else "_nolegend_"
label2 = "True future values" if j==0 else "_nolegend_"
label3 = "Predictions" if j==0 else "_nolegend_"
label4 = "Ruitu values" if j==0 else "_nolegend_"
plt.plot(range(len(past)), past, "o-g",
label=label1)
plt.plot(range(len(past),
len(true)+len(past)), true, "x--g", label=label2)
plt.plot(range(len(past), len(pred)+len(past)), pred, "o--y",
label=label3)
if input_ruitu is not None:
plt.plot(range(len(past), len(ruitu)+len(past)), ruitu, "o--r",
label=label4)
plt.legend(loc='best')
plt.title("Predictions v.s. true values v.s. Ruitu")
plt.show()
class RNN_Class(WeatherConv1D):
def __init__(self, num_output_features, num_decoder_features,
target_sequence_length,
num_steps_to_predict, regulariser = None,
lr=0.001, decay=0, loss = "mse",
layers=[35, 35]):
self.num_output_features = num_output_features
self.num_decoder_features = num_decoder_features
self.target_sequence_length = target_sequence_length
self.num_steps_to_predict = num_steps_to_predict
self.regulariser = regulariser
self.layers = layers
self.lr = lr
self.decay = decay
self.loss = loss
self.pred_result = None
#self.batch_size = batch_size
print('Initialized!')
def build_graph(self):
keras.backend.clear_session() # clear session/graph
self.model = RNN_builder(self.num_output_features, self.num_decoder_features,
self.target_sequence_length,
self.num_steps_to_predict, self.regulariser,
self.lr, self.decay, self.loss, self.layers)
print(self.model.summary())
class Seq2Seq_Class(Enc_Dec):
def __init__(self, id_embd, time_embd,
num_input_features, num_output_features, num_decoder_features,
input_sequence_length, target_sequence_length,
num_steps_to_predict, regulariser = None,
lr=0.001, decay=0, loss = "mse",
layers=[35, 35], model_save_path='../models',
model_structure_name='seq2seq_model.json', model_weights_name='seq2seq_model_weights.h5'):
super().__init__(num_input_features, num_output_features, num_decoder_features,
input_sequence_length, target_sequence_length,
num_steps_to_predict, regulariser = None,
lr=lr, decay=decay, loss = loss,
layers=layers)
self.id_embd = id_embd
self.time_embd = time_embd
self.val_loss_list=[]
self.train_loss_list=[]
self.current_mean_val_loss = None
self.early_stop_limit = 10 # with the unit of Iteration Display
self.EARLY_STOP=False
self.pred_var_result = []
self.pi_dic={0.95:1.96, 0.9:1.645, 0.8:1.28, 0.68:1.}
self.target_list=['t2m','rh2m','w10m']
self.obs_range_dic={'t2m':[-30,42], # Official value: [-20,42]
'rh2m':[0.0,100.0],
'w10m':[0.0, 30.0]}
self.obs_and_output_feature_index_map = {'t2m':0,'rh2m':1,'w10m':2}
self.ruitu_feature_index_map = {'t2m':1,'rh2m':3,'w10m':4}
self.model_save_path = model_save_path
self.model_structure_name=model_structure_name
self.model_weights_name=model_weights_name
def build_graph(self):
#keras.backend.clear_session() # clear session/graph
self.optimizer = keras.optimizers.Adam(lr=self.lr, decay=self.decay)
self.model = Seq2Seq_MVE_subnets_swish(id_embd=True, time_embd=True,
lr=self.lr, decay=self.decay,
num_input_features=self.num_input_features, num_output_features=self.num_output_features,
num_decoder_features=self.num_decoder_features, layers=self.layers,
loss=self.loss, regulariser=self.regulariser)
def _mve_loss(y_true, y_pred):
pred_u = crop(2,0,3)(y_pred)
pred_sig = crop(2,3,6)(y_pred)
print(pred_sig)
#exp_sig = tf.exp(pred_sig) # avoid pred_sig is too small such as zero
#precision = 1./exp_sig
precision = 1./pred_sig
#log_loss= 0.5*tf.log(exp_sig)+0.5*precision*((pred_u-y_true)**2)
log_loss= 0.5*tf.log(pred_sig)+0.5*precision*((pred_u-y_true)**2)
log_loss=tf.reduce_mean(log_loss)
return log_loss
print(self.model.summary())
self.model.compile(optimizer = self.optimizer, loss=_mve_loss)
def sample_batch(self, data_inputs, ground_truth, ruitu_inputs, batch_size, certain_id=None, certain_feature=None):
max_i, _, max_j, _ = data_inputs.shape # Example: (1148, 37, 10, 9)-(sample_ind, timestep, sta_id, features)
id_ = np.random.randint(max_j, size=batch_size)
i = np.random.randint(max_i, size=batch_size)
batch_inputs = data_inputs[i,:,id_,:]
batch_ouputs = ground_truth[i,:,id_,:]
batch_ruitu = ruitu_inputs[i,:,id_,:]
# id used for embedding
if self.id_embd and (not self.time_embd):
expd_id = np.expand_dims(id_,axis=1)
batch_ids = np.tile(expd_id,(1,37))
return batch_inputs, batch_ruitu, batch_ouputs, batch_ids
elif (not self.id_embd) and (self.time_embd):
time_range = np.array(range(37))
batch_time = np.tile(time_range,(batch_size,1))
#batch_time = np.expand_dims(batch_time, axis=-1)
return batch_inputs, batch_ruitu, batch_ouputs, batch_time
elif (self.id_embd) and (self.time_embd):
expd_id = np.expand_dims(id_,axis=1)
batch_ids = np.tile(expd_id,(1,37))
time_range = np.array(range(37))
batch_time = np.tile(time_range,(batch_size,1))
#batch_time = np.expand_dims(batch_time, axis=-1)
return batch_inputs, batch_ruitu, batch_ouputs, batch_ids, batch_time
elif (not self.id_embd) and (not self.time_embd):
return batch_inputs, batch_ruitu, batch_ouputs
def fit(self, train_input_obs, train_input_ruitu, train_labels,
val_input_obs, val_input_ruitu, val_labels, val_ids, val_times, batch_size,
iterations=300, validation=True):
print('Train batch size: {}'.format(batch_size))
print('Validation on data size of {};'.format(val_input_obs.shape[0]))
early_stop_count = 0
for i in range(iterations):
batch_inputs, batch_ruitu, batch_labels, batch_ids, batch_time = self.sample_batch(train_input_obs, train_labels,
train_input_ruitu, batch_size=batch_size)
#batch_placeholders = np.zeros_like(batch_labels)
loss_ = self.model.train_on_batch(x=[batch_inputs, batch_ruitu, batch_ids, batch_time],
y=[batch_labels])
if (i+1)%50 == 0:
print('Iteration:{}/{}. Training batch MLE loss:{}'.
format(i+1, iterations, loss_))
if validation :
self.evaluate(val_input_obs, val_input_ruitu, val_labels, val_ids, val_times, each_station_display=False)
if len(self.val_loss_list) >0: # Early stopping
if(self.current_mean_val_loss) <= min(self.val_loss_list): # compare with the last early_stop_limit values except SELF
early_stop_count = 0
model_json = self.model.to_json()
with open(self.model_save_path+self.model_structure_name, "w") as json_file:
json_file.write(model_json)
self.model.save_weights(self.model_save_path+self.model_weights_name)
else:
early_stop_count +=1
print('Early-stop counter:', early_stop_count)
if early_stop_count == self.early_stop_limit:
self.EARLY_STOP=True
break
print('###'*10)
if self.EARLY_STOP:
print('Loading the best model before early-stop ...')
self.model.load_weights(self.model_save_path+self.model_weights_name)
print('Training finished! Detailed val MLE loss:')
self.evaluate(val_input_obs, val_input_ruitu, val_labels, val_ids, val_times, each_station_display=True)
def evaluate(self, data_input_obs, data_input_ruitu, data_labels, data_ids, data_time, each_station_display=False):
all_loss=[]
for i in range(10): # iterate for each station. (sample_ind, timestep, staionID, features)
#batch_placeholders = np.zeros_like(data_labels[:,:,i,:])
val_loss= self.model.evaluate(x=[data_input_obs[:,:,i,:], data_input_ruitu[:,:,i,:], data_ids[:,:,i], data_time],
y=[data_labels[:,:,i,:]], verbose=False)
all_loss.append(val_loss)
if each_station_display:
print('\tFor station 9000{}, val MLE loss: {}'.format(i+1, val_loss))
self.current_mean_val_loss = np.mean(all_loss)
print('Mean val MLE loss:', self.current_mean_val_loss)
self.val_loss_list.append(self.current_mean_val_loss)
def predict(self, batch_inputs, batch_ruitu, batch_ids, batch_times):
'''
Input:
Output:
pred_result (mean value) : (None, 10,37,3). i.e., (sample_nums, stationID, timestep, features)
pred_var_result (var value) : (None, 10,37,3)
'''
pred_result_list = []
pred_var_list = []
#pred_std_list =[]
for i in range(10):
result = self.model.predict(x=[batch_inputs[:,:,i,:], batch_ruitu[:,:,i,:], batch_ids[:,:,i], batch_times])
var_result = result[:,:,3:6] # Variance
result = result[:,:,0:3] # Mean
#result = np.squeeze(result, axis=0)
pred_result_list.append(result)
#var_result = np.squeeze(var_result, axis=0)
pred_var_list.append(var_result)
pred_result = np.stack(pred_result_list, axis=1)
pred_var_result = np.stack(pred_var_list, axis=1)
print('Predictive shape (None, 10,37,3) means (sample_nums, stationID, timestep, features). \
Features include: t2m, rh2m and w10m')
self.pred_result = pred_result
self.pred_var_result = pred_var_result
#self.pred_std_result = np.sqrt(np.exp(self.pred_var_result[:,:,i,j])) # Calculate standard deviation
return pred_result, pred_var_result
def renorm_for_visualization(self, obs_inputs, ruitu_inputs, pred_mean_result, pred_var_result, ground_truth=None):
'''
obs_inputs: (None, 28, 10, 9)
ruitu_inputs: (None, 37, 10, 29)
pred_mean_result: (None, 10, 37, 3)
pred_var_result: (None, 10, 37, 3)
ground_truth: (None, 37, 10, 3)
#self.target_list=['t2m','rh2m','w10m']
#self.obs_range_dic={'t2m':[-30,42],
# 'rh2m':[0.0,100.0],
# 'w10m':[0.0, 30.0]}
#self.obs_and_output_feature_index_map = {'t2m':0,'rh2m':1,'w10m':2}
#self.ruitu_feature_index_map = {'t2m':1,'rh2m':3,'w10m':4}
#TODO:
'''
for target_v in self.target_list:
temp1 = obs_inputs[:,:,:,self.obs_and_output_feature_index_map[target_v]]
temp2 = ruitu_inputs[:,:,:,self.ruitu_feature_index_map[target_v]]
temp3 = pred_mean_result[:,:,:,self.obs_and_output_feature_index_map[target_v]]
#temp4 = pred_var_result[:,:,:,self.obs_and_output_feature_index_map[target_v]]
obs_inputs[:,:,:,self.obs_and_output_feature_index_map[target_v]] = renorm(temp1, self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
ruitu_inputs[:,:,:,self.ruitu_feature_index_map[target_v]] = renorm(temp2, self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
pred_mean_result[:,:,:,self.obs_and_output_feature_index_map[target_v]] = renorm(temp3, self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
#pred_var_result[:,:,:,self.obs_and_output_feature_index_map[target_v]] = renorm(temp4, self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
if ground_truth is not None:
temp5 = ground_truth[:,:,:,self.obs_and_output_feature_index_map[target_v]]
ground_truth[:,:,:,self.obs_and_output_feature_index_map[target_v]] = renorm(temp5, self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
if ground_truth is not None:
return obs_inputs, ruitu_inputs, pred_mean_result, pred_var_result, ground_truth
else:
return obs_inputs, ruitu_inputs, pred_mean_result, pred_var_result
def calc_uncertainty_info(self, verbose=False):
'''
Verbose: Display uncertainty for each feature i.e., (t2m, rh2m, w10m)
#TODO: Refactor the double 'for' part.
'''
assert len(self.pred_var_result)>0, 'Error! You must run predict() before running calc_uncertainty_info()'
print('The uncertainty info are calculated on {} predicted samples with shape {}'
.format(len(self.pred_var_result), self.pred_var_result.shape))
#
if verbose:
assert self.target_list == ['t2m','rh2m','w10m'], 'ERROR, list changed!'
for j, target_v in enumerate(['t2m','rh2m','w10m']):
print('For feature {}:'.format(target_v))
for i in range(37):
#unctt_var = np.exp(self.pred_var_result[:,:,i,j])
unctt_std = np.sqrt(unctt_var)
unctt_mean_std = np.mean(unctt_std)
unctt_mean_var = np.mean(unctt_var)
#renorm_unctt_mean_std = renorm(unctt_mean_std, self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
print('\tTime:{}-Variance:{:.4f}; Std:{:.4f};'.
format(i+1, unctt_mean_var, unctt_mean_std))
else:
for i in range(37):
unctt_var = np.exp(self.pred_var_result[:,:,i,:])
unctt_std = np.sqrt(unctt_var)
unctt_mean_std = np.mean(unctt_std)
unctt_mean_var = np.mean(unctt_var)
#renorm_unctt_mean_std = 0
print('Time:{}-Variance:{:.4f}; Std:{:.4f};'.
format(i+1, unctt_mean_var, unctt_mean_std))
def minus_plus_std_strategy(self, pred_mean, pred_var, feature_name,\
timestep_to_ensemble=21, alpha=0):
'''
This stratergy aims to calculate linear weighted at specific timestep (timestep_to_ensemble) between prediction and ruitu as formula:
(alpha)*pred_mean + (1-alpha)*ruitu_inputs
pred_mean: (10, 37, 3)
pred_var: (10, 37, 3)
timestep_to_ensemble: int32 (From 0 to 36)
'''
print('Using minus_plus_var_strategy with alpha {}'.format(alpha))
assert 0<=timestep_to_ensemble<=36 , 'Please ensure 0<=timestep_to_ensemble<=36!'
assert -0.3<= alpha <=0.3, '-0.3<= alpha <=0.3!'
assert pred_mean.shape == (10, 37, 3), 'Error! This funtion ONLY works for \
one data sample with shape (10, 37, 3). Any data shape (None, 10, 37, 3) will leads this error!'
pred_std = np.sqrt(np.exp(pred_var))
print('alpha:',alpha)
pred_mean[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]] = \
pred_mean[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]] + \
alpha * pred_std[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]]
return pred_mean
def linear_ensemble_strategy(self, pred_mean, pred_var, ruitu_inputs, feature_name,\
timestep_to_ensemble=21, alpha=1):
'''
This stratergy aims to calculate linear weighted at specific timestep (timestep_to_ensemble) between prediction and ruitu as formula:
(alpha)*pred_mean + (1-alpha)*ruitu_inputs
pred_mean: (10, 37, 3)
pred_var: (10, 37, 3)
ruitu_inputs: (37,10,29). Need Swamp to(10,37,29) FIRSTLY!!
timestep_to_ensemble: int32 (From 0 to 36)
'''
assert 0<= alpha <=1, 'Please ensure 0<= alpha <=1 !'
assert pred_mean.shape == (10, 37, 3), 'Error! This funtion ONLY works for \
one data sample with shape (10, 37, 3). Any data shape (None, 10, 37, 3) will leads this error!'
#pred_std = np.sqrt(np.exp(pred_var))
ruitu_inputs = np.swapaxes(ruitu_inputs,0,1)
print('alpha:',alpha)
pred_mean[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]] = \
(alpha)*pred_mean[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]] + \
(1-alpha)*ruitu_inputs[:,timestep_to_ensemble:, self.ruitu_feature_index_map[feature_name]]
print('Corrected pred_mean shape:', pred_mean.shape)
return pred_mean
def fuzzy_ensemble_strategy(self, pred_mean, pred_var, feature_name,\
timestep_to_ensemble=21, alpha=0):
'''
This stratergy aims to calculate linear weighted at specific timestep (timestep_to_ensemble) between prediction and ruitu as formula:
(alpha)*pred_mean + (1-alpha)*ruitu_inputs
pred_mean: (10, 37, 3)
pred_var: (10, 37, 3)
timestep_to_ensemble: int32 (From 0 to 36)
'''
print('Using fuzzy_ensemble_strategy with alpha {}'.format(alpha))
assert 0<=timestep_to_ensemble<=36 , 'Please ensure 0<=timestep_to_ensemble<=36!'
assert -0.4<= alpha <=0.4, 'Please ensure -0.4<= alpha <=0.4 !'
assert pred_mean.shape == (10, 37, 3), 'Error! This funtion ONLY works for \
one data sample with shape (10, 37, 3). Any data shape (None, 10, 37, 3) will leads this error!'
pred_std = np.sqrt(np.exp(pred_var))
#print('normalizing for Std. after timestep:', timestep_to_ensemble)
temp_std = pred_std[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]]
norm_std = temp_std / np.max(temp_std)
#print('norm_std shape', norm_std.shape)
dim_0, dim_1 = norm_std.shape
reshaped_std = norm_std.reshape(-1)
from skfuzzy import trimf
fuzzy_degree = trimf(reshaped_std, [0., 1, 1.2])
fuzzy_degree = fuzzy_degree.reshape(dim_0, dim_1)
#print('fuzzy_degree shape:',fuzzy_degree.shape)
#print('temp_std shape:',temp_std.shape)
pred_mean[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]] = \
pred_mean[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]] + \
fuzzy_degree*alpha*temp_std
#pred_mean[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]] + \
#alpha * pred_std[:,timestep_to_ensemble:,self.obs_and_output_feature_index_map[feature_name]]
#print('pred_mean.shape',pred_mean.shape)
return pred_mean
pass
def renorm_for_submit(self, pred_mean, pred_var, ruitu_inputs, timestep_to_ensemble=21, alpha=1):
'''
Overwrite for Seq2Seq_MVE Class
pred_mean: shape of (10, 37, 3)
pred_var: shape of (10, 37, 3)
ruitu_inputs: shape of (10, 37, 3)
timestep_to_ensemble: int32 (From 0 to 36)
# TODO: Add three strategies for output
'''
assert self.pred_result is not None, 'You must run self.predict(batch_inputs, batch_ruitu) firstly!!'
assert pred_mean.shape == (10, 37, 3), 'Error! This funtion ONLY works for one data sample with shape (10, 37, 3). Any data shape (None, 10, 37, 3) will leads this error!'
df_empty = pd.DataFrame(columns=['FORE_data', 't2m', 'rh2m', 'w10m'])
for j, target_v in enumerate(self.target_list):
series_ids = pd.Series()
series_targets = pd.Series()
#print('Feature {}, timestep_to_ensemble: {}, weighted alpha: {}'.
# format(target_v, timestep_to_ensemble, alpha))
#pred_mean = self.linear_ensemble_strategy(pred_mean, pred_var,
# ruitu_inputs, target_v, timestep_to_ensemble, alpha)
#pred_mean =self.minus_plus_std_strategy(pred_mean, pred_var, target_v,\
# timestep_to_ensemble, alpha)
#pred_mean = self.fuzzy_ensemble_strategy(pred_mean, pred_var, target_v,\
# timestep_to_ensemble, alpha=0.)
renorm_value = renorm(pred_mean[:,:,j], self.obs_range_dic[target_v][0], self.obs_range_dic[target_v][1])
for i in range(10):
if i != 9:
id_num = '0'+str(i+1)
else:
id_num = str(10)
sta_name_time = '900'+id_num+'_'
time_str_list=[]
for t in range(37):
if t < 10:
time_str= sta_name_time + '0'+ str(t)
else:
time_str = sta_name_time + str(t)
time_str_list.append(time_str)
series_id = pd.Series(time_str_list)
series_target = pd.Series(renorm_value[i])
series_ids = pd.concat([series_ids, series_id])
series_targets = pd.concat([series_targets, series_target])
df_empty['FORE_data'] = series_ids
df_empty[target_v] = series_targets
return df_empty
def plot_prediction(self, x, y_true, y_pred, intervals=None, input_ruitu=None, pi_degree=0.8, renorm_flag=False):
"""Plots the predictions.
Arguments
---------
x: Input sequence of shape (input_sequence_length,
dimension_of_signal) E.g. (28, 1)
y_true: True output sequence of shape (input_sequence_length,
dimension_of_signal) E.g. (35, 1)
y_pred: Predicted output sequence (input_sequence_length,
dimension_of_signal) E.g. (35, 1)
input_ruitu: Ruitu output sequence E.g. (35, 1)
pi_degree: Confident Level such as 0.95, 0.9, 0.8, and 0.68 etc.
"""
plt.figure(figsize=(12, 3))
output_dim = x.shape[-1]# feature dimension
for j in range(output_dim):
past = x[:, j]
true = y_true[:, j]
pred = y_pred[:, j]
if input_ruitu is not None:
ruitu = input_ruitu[:, j]
if intervals is not None:
pi_var = intervals[:, j]
pi_var = np.sqrt(np.exp(pi_var))
label1 = "Seen (past) values" if j==0 else "_nolegend_"
label2 = "True future values" if j==0 else "_nolegend_"
label3 = "Predictions" if j==0 else "_nolegend_"
label4 = "Ruitu values" if j==0 else "_nolegend_"
label5 = "Lower-Upper bound" if j==0 else "_nolegend_"
plt.plot(range(len(past)), past, "o-g",
label=label1)
plt.plot(range(len(past),
len(true)+len(past)), true, "x--g", label=label2)
plt.plot(range(len(past), len(pred)+len(past)), pred, ".--b",
label=label3)
if input_ruitu is not None:
plt.plot(range(len(past), len(ruitu)+len(past)), ruitu, ".--r",
label=label4)
if intervals is not None:
#print(intervals.shape)
print(pi_var.shape)
up_bound = pred + self.pi_dic[pi_degree]*pi_var
low_bound = pred - self.pi_dic[pi_degree]*pi_var
plt.fill_between(range(len(past), len(ruitu)+len(past)),
up_bound, low_bound, facecolor='blue', alpha=0.1)
plt.legend(loc='best')
plt.title("Predictions v.s. true values v.s. Ruitu")
plt.show()
class Enc_Dec_Embd(Enc_Dec):
def build_graph(self):
keras.backend.clear_session() # clear session/graph
self.optimiser = keras.optimizers.Adam(lr=self.lr, decay=self.decay)
# Define an input sequence.
encoder_inputs = keras.layers.Input(shape=(None, self.num_input_features), name='encoder_inputs')
# Create a list of RNN Cells, these are then concatenated into a single layer
# with the RNN layer.
encoder_cells = []
for hidden_neurons in self.layers:
encoder_cells.append(keras.layers.GRUCell(hidden_neurons,
kernel_regularizer = self.regulariser,
recurrent_regularizer = self.regulariser,
bias_regularizer = self.regulariser))
encoder = keras.layers.RNN(encoder_cells, return_state=True)
encoder_outputs_and_states = encoder(encoder_inputs)
# Discard encoder outputs and only keep the states.
encoder_states = encoder_outputs_and_states[1:]
# Define a decoder sequence.
decoder_inputs = keras.layers.Input(shape=(None, self.num_decoder_features), name='decoder_inputs')
decoder_inputs_id = keras.layers.Input(shape=(None,), name='id_inputs')
decoder_inputs_id_embd = Embedding(input_dim=10, output_dim=2, name='id_embedding')(decoder_inputs_id)
#decoder_inputs_time = keras.layers.Input(shape=(None,), name='time_inputs')
#decoder_inputs_time_embd = Embedding(input_dim=37, output_dim=2, name='time_embedding')(decoder_inputs_time)
decoder_concat = concatenate([decoder_inputs, decoder_inputs_id_embd], axis=-1)
decoder_cells = []
for hidden_neurons in self.layers:
decoder_cells.append(keras.layers.GRUCell(hidden_neurons,
kernel_regularizer = self.regulariser,
recurrent_regularizer = self.regulariser,
bias_regularizer = self.regulariser))
decoder = keras.layers.RNN(decoder_cells, return_sequences=True, return_state=True)
decoder_outputs_and_states = decoder(decoder_concat, initial_state=encoder_states)
decoder_outputs = decoder_outputs_and_states[0]
#decoder_dense1 = keras.layers.Dense(units=32,
# activation='relu',
# kernel_regularizer = self.regulariser,
# bias_regularizer = self.regulariser, name='dense_relu')
output_dense = keras.layers.Dense(self.num_output_features,
activation='sigmoid',
kernel_regularizer = self.regulariser,
bias_regularizer = self.regulariser, name='output_sig')
#densen1=decoder_dense1(decoder_outputs)
decoder_outputs = output_dense(decoder_outputs)
# Create a model using the functional API provided by Keras.
self.model = keras.models.Model(inputs=[encoder_inputs, decoder_inputs, decoder_inputs_id], outputs=decoder_outputs)
self.model.compile(optimizer = self.optimiser, loss=self.loss)
print(self.model.summary())
def sample_batch(self, data_inputs, ground_truth, ruitu_inputs, batch_size, certain_id=None, certain_feature=None):
max_i, _, max_j, _ = data_inputs.shape # Example: (1148, 37, 10, 9)-(sample_ind, timestep, sta_id, features)
if certain_id == None and certain_feature == None:
id_ = np.random.randint(max_j, size=batch_size)
i = np.random.randint(max_i, size=batch_size)
batch_inputs = data_inputs[i,:,id_,:]
batch_ouputs = ground_truth[i,:,id_,:]
batch_ruitu = ruitu_inputs[i,:,id_,:]
# id used for embedding
expd_id = np.expand_dims(id_,axis=1)
batch_ids = np.tile(expd_id,(1,37))
#batch_time =
elif certain_id != None:
pass
return batch_inputs, batch_ruitu, batch_ouputs, batch_ids
def fit(self, train_input_obs, train_input_ruitu, train_labels,
val_input_obs, val_input_ruitu, val_labels, val_ids, batch_size,
iterations=300, validation=True):
print('Train batch size: {}'.format(batch_size))
print('Validation on data size of {};'.format(val_input_obs.shape[0]))
for i in range(iterations):
batch_inputs, batch_ruitu, batch_labels, batch_ids = self.sample_batch(train_input_obs, train_labels,
train_input_ruitu, batch_size=batch_size)
#print(batch_ids.shape)
loss_ = self.model.train_on_batch(x=[batch_inputs, batch_ruitu, batch_ids],
y=[batch_labels])
if (i+1)%50 == 0:
print('Iteration:{}/{}. Training batch loss:{}'.
format(i+1, iterations, loss_))
if validation :
self.evaluate(val_input_obs, val_input_ruitu, val_labels, val_ids, each_station_display=False)
print('###'*10)
print('Train finish! Total validation loss:')
self.evaluate(val_input_obs, val_input_ruitu, val_labels, val_ids, each_station_display=True)
def evaluate(self, data_input_obs, data_input_ruitu, data_labels, data_ids, each_station_display=False, each_feature_display=False):
#assert data_input_ruitu.shape[0] == data_input_obs.shape[0] == data_labels.shape[0], 'Shape Error'
#assert data_input_obs.shape[1] == 28 and data_input_obs.shape[2] == 10 and data_input_obs.shape[3] == 9, 'Error! Obs input shape must be (None, 28,10,9)'
#assert data_input_ruitu.shape[1] == 37 and data_input_ruitu.shape[2] == 10 and data_input_ruitu.shape[3] == 29, 'Error! Ruitu input shape must be (None, 37,10,29)'
#assert data_labels.shape[1] == 37 and data_labels.shape[2] == 10 and data_labels.shape[3] == 3, 'Error! Ruitu input shape must be (None, 37,10,3)'
all_loss=[]
for i in range(10): # iterate for each station. (sample_ind, timestep, staionID, features)
#print(data_ids.shape)
val_loss= self.model.evaluate(x=[data_input_obs[:,:,i,:], data_input_ruitu[:,:,i,:], data_ids[:,:,i]],
y=data_labels[:,:,i,:], verbose=False)
all_loss.append(val_loss)
if each_station_display:
print('\tFor station 9000{}, evaluated loss: {}'.format(i+1, val_loss))
print('Mean evaluated loss on all stations:', np.mean(all_loss))
def predict(self, batch_inputs, batch_ruitu, batch_ids):
#assert batch_ruitu.shape[0] == batch_inputs.shape[0], 'Shape Error'
#assert batch_inputs.shape[1] == 28 and batch_inputs.shape[2] == 10 and batch_inputs.shape[3] == 9, 'Error! Obs input shape must be (None, 28,10,9)'
#assert batch_ruitu.shape[1] == 37 and batch_ruitu.shape[2] == 10 and batch_ruitu.shape[3] == 29, 'Error! Ruitu input shape must be (None, 37,10, 29)'
pred_result_list = []
for i in range(10):
#print('Predict for station: 9000{}'.format(i+1))
result = self.model.predict(x=[batch_inputs[:,:,i,:], batch_ruitu[:,:,i,:], batch_ids[:,:,i]])
result = np.squeeze(result, axis=0)
#all_pred[i] = result
pred_result_list.append(result)
#pass
pred_result = np.stack(pred_result_list, axis=0)
#return all_pred, pred_result
print('Predict shape (10,37,3) means (stationID, timestep, features). Features include: t2m, rh2m and w10m')
self.pred_result = pred_result
return pred_result
def renorm(norm_value, min_v ,max_v):
real_v = norm_value * (max_v-min_v) + min_v
return real_v
