<jupyter_start><jupyter_text>Previous code for data preprocessing<jupyter_code>gpu_info = !nvidia-smi
gpu_info = '\n'.join(gpu_info)
if gpu_info.find('failed') >= 0:
  print('Select the Runtime → "Change runtime type" menu to enable a GPU accelerator, ')
  print('and then re-execute this cell.')
else:
  print(gpu_info)
# Setting up the environment
! pip install numpy scipy matplotlib xarray pandas netcdf4 tqdm tensorflow scikit-learn seaborn
import s3fs
import seaborn as sns
import pandas as pd
import pyarrow
import re
import random
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from glob import glob
from os.path import join, exists
from matplotlib.colors import LogNorm
import tensorflow as tf
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from tensorflow.keras.layers import Dense, Input
from tensorflow.keras.models import Model
from sklearn.model_selection import GroupKFold
from tensorflow.keras.optimizers import Adam
%matplotlib inline

seed = 8886
np.random.seed(seed)
random.seed(seed)
tf.random.set_seed(seed)
# Load functions. Copied from the instruction notebook.
# -----------------------------------------------
""" Data Processing """

def fetch_data(f):
    """ Load data directly from Amazon S3 storage """
    
    obj = fs.open(f)
    df = pd.read_parquet(obj)
    return df

# def prepare_data(data, input_columns, output_columns):
#     """ Splits data to input/output with corresponding lag """
    
#     data['id'] = data['id'].apply(lambda x: int(x[3:])) # change expnumber to int for sorting
#     data.index.name = 'indx'
#     data = data.sort_values(['id','indx'])
#     x_data = data[input_columns].iloc[:-1, :]
#     y_data = data[output_columns].iloc[1:, :]
    
#     return x_data, y_data

def prepare_data(data, input_columns, output_columns):
    """ Splits data to input/output with corresponding lag """
    data['id'] = data['id'].apply(lambda x: int(x[3:])) # change expnumber to int for sorting
    data.index.name = 'indx'
    data = data.sort_values(['id','indx'])
    # lag in/out by experiment
    mask = data['Time [s]'] != data['Time [s]'].max()
    x_data = data.loc[mask, input_columns]
    mask = data['Time [s]'] != data['Time [s]'].min()
    y_data = data.loc[mask, output_columns]
    return x_data, y_data

def get_starting_conds(data, input_vars, exp, starting_ts=0):
    """ Take data and expiriment number to gather initial starting condition for box emulator"""
    exp_data = data[data['id'] == exp]
    
    if exp_data.shape[0] == 0:
        raise ValueError('That expiriment number does not exist in this dataset.')
    
    ts_data = exp_data.iloc[starting_ts:starting_ts+1, :]

    return ts_data
# -----------------------------------------------
""" Evaluation """

def calc_pdf_hist(x, x_bins):
    return np.histogram(x, x_bins, density=True)[0]

def rmse(y_true, y_pred):
    return np.sqrt(mean_squared_error(y_true, y_pred))

def hellinger(x, pdf_p, pdf_q):
    pdf_distances = (np.sqrt(pdf_p) - np.sqrt(pdf_q)) ** 2
    return np.trapz(pdf_distances, x) / 2

def hellinger_distance(y_true, y_pred, bins=50):
    bin_points = np.linspace(np.minimum(y_true.min(), y_pred.min()),
                       np.maximum(y_true.max(), y_pred.max()),
                       bins)
    bin_centers = 0.5 * (bin_points[:-1] + bin_points[1:])
    y_true_pdf = calc_pdf_hist(y_true, bin_points)
    y_pred_pdf = calc_pdf_hist(y_pred, bin_points)
    return hellinger(bin_centers, y_true_pdf, y_pred_pdf)

def r2_corr(y_true, y_pred):
    return np.corrcoef(y_true, y_pred)[0, 1] ** 2

def evaluate_mod(true_output, model_results):
    
    print('RMSE: Precursor: {:.5f}, Gas: {:.5f}, Aerosols: {:.5f}'.format(
        rmse(true_output.iloc[:,0], model_results.iloc[:,0]),
        rmse(true_output.iloc[:,1], model_results.iloc[:,1]),
        rmse(true_output.iloc[:,2], model_results.iloc[:,2])))
    print('R2: Precursor: {:.5f}, Gas: {:.5f}, Aerosols: {:.5f}'.format(
        r2_corr(true_output.iloc[:,0], model_results.iloc[:,0]),
        r2_corr(true_output.iloc[:,1], model_results.iloc[:,1]),
        r2_corr(true_output.iloc[:,2], model_results.iloc[:,2])))
    print('Hellenger Distance: Precursor: {:.5f}, Gas: {:.5f}, Aerosols: {:.5f}'.format(
        hellinger_distance(true_output.iloc[:,0], model_results.iloc[:,0]),
        hellinger_distance(true_output.iloc[:,1], model_results.iloc[:,1]),
        hellinger_distance(true_output.iloc[:,2], model_results.iloc[:,2])))
    
    return
# -----------------------------------------------
""" Models """
def dense_neural_net(num_layers=2, num_neurons=100, activation="relu", learning_rate=0.001,
                     num_input_columns=9, num_output_columns=3):
    """ Build and return dense neural network with specified hyperparameters"""
    
    input_layer = Input(num_input_columns)
    n_net = input_layer
    for l in range(num_layers):
        n_net= Dense(num_neurons, activation=activation)(n_net)
    out = Dense(num_output_columns)(n_net)
    mod = Model(input_layer, out)
    mod.compile(Adam(learning_rate=learning_rate), "mse")
    
    return mod

def box_emulator(neural_network, starting_conds, input_scaler, output_scaler, num_timesteps=1439):
    """ Forward box emulator
    Args:
        neural_netwrok: NN model object that predits output at t+1
        starting_conds: Initial input conditions used for first prediction
        input_scaler: Input scaler object
        onput_scaler: Onput scaler object
        num_timesteps: how many timesteps forward to run emulator
    Returns:
        results: Pandas DataFrame of box emulator output
    """
    mod = neural_network
    scaled_input = input_scaler.transform(starting_conds.iloc[:,1:-1])
    static_input = scaled_input[:,-6:]

    for i in range(num_timesteps):

        if i == 0:

            pred = mod.predict(scaled_input)
            new_input = np.concatenate([pred,static_input], axis=1)
            pred_array = pred

        else:

            pred = mod.predict(new_input)
            new_input = np.concatenate([pred,static_input], axis=1)
            pred_array = np.concatenate([pred_array, pred], axis=0)
            
    results = pd.DataFrame(output_scaler.inverse_transform(pred_array))
    results['id'] = starting_conds.iloc[0,-1]
    results.columns = starting_conds.columns[[1,2,3,-1]]
    return results

# def run_emulator_ensemble(input_data, mod, num_expiriments, input_vars, input_scaler, 
#                           output_scaler, run_length=1439):
#     """ Run an ensemble of emulators on n amount of random expiriments
#     Args:
#         input_data: unscaled dataframe collection of expiriment data
#         mod: NN model object that predits output at t+1
#         num_expiriments: How many expiriments to run
#         input_vars: Input variables for model
#         output_vars: Output variables for model
#         run_length: number of timesteps to run each emulator forward
#     Returns:
#         pred_df: Aggregated Pandas DataFrame of predictions for each expiriment
#         true_df: Corresponding Dataframe with raw expiriment data
#     """
    output_cols = input_data.columns[[1,2,3, -1]]
    exps = sorted(random.sample(list(input_data['id'].unique()), num_expiriments))
    pred_df  = pd.DataFrame(columns = output_cols)
    true_df = input_data.loc[input_data['id'].isin(exps)][output_cols].reset_index(drop=True)
    true_df =  true_df.groupby('id').apply(lambda x: x.iloc[1:,:]).reset_index(drop=True)  
    
    for i in range(num_expiriments):
        
        init_conds = get_starting_conds(input_data, input_vars, exps[i])
        preds = box_emulator(mod, init_conds, input_scaler, output_scaler, run_length)
        preds['id'] = exps[i]
        preds.columns = pred_df.columns
        pred_df = pred_df.append(preds, ignore_index=True)

    return pred_df, true_df

def run_emulator_ensemble(input_data, output_data, mod, num_expiriments, input_vars, input_scaler, 
                          output_scaler, run_length=1439):
    """ Run an ensemble of emulators on n amount of random expiriments
    Args:
        input_data: unscaled dataframe collection of expiriment data
        mod: NN model object that predits output at t+1
        num_expiriments: How many expiriments to run
        input_vars: Input variables for model
        output_vars: Output variables for model
        run_length: number of timesteps to run each emulator forward
    Returns:
        pred_df: Aggregated Pandas DataFrame of predictions for each expiriment
        true_df: Corresponding Dataframe with raw expiriment data
    """
    output_cols = input_data.columns[[1,2,3, -1]]
    exps = sorted(random.sample(list(input_data['id'].unique()), num_expiriments))
    pred_df  = pd.DataFrame(columns = output_cols)
    true_df = output_data.loc[output_data['id'].isin(exps)].iloc[:,1:-1].reset_index(drop=True)
     
    for i in range(num_expiriments):
        init_conds = get_starting_conds(input_data, input_vars, exps[i])
        preds = box_emulator(mod, init_conds, input_scaler, output_scaler, run_length)
        preds['id'] = exps[i]
        preds.columns = pred_df.columns
        pred_df = pred_df.append(preds, ignore_index=True)
    return pred_df, true_df<jupyter_output><empty_output><jupyter_text>Load the data<jupyter_code># load data
fs = s3fs.S3FileSystem(anon=True)
gecko_files = fs.ls("ncar-aiml-data-commons/gecko/")[1:]
test, train, val = map(fetch_data, gecko_files)
train.head()
# get some summary statistics of the training data
train.groupby('id').mean().describe()
# Data Preperation
'''
Copied from the instruction notebook

Here we prepare the data for machine learning by taking the respective variables 
from each column, offsetting the output data by 1 timestep (this is done inside 
the `prepare_data()` function, and scale. Make sure to not re-fit the scaler on 
the validation/test data and only `transform()` it. 
'''

input_vars, output_vars = train.columns, train.columns[[0,1,2,3,-1]]
train_in, train_out = prepare_data(train, input_vars, output_vars)
val_in, val_out = prepare_data(val, input_vars, output_vars)<jupyter_output><empty_output><jupyter_text>Data scaling<jupyter_code># Min max scaling
input_scaler = MinMaxScaler()
output_scaler = MinMaxScaler()

scaled_train_in = input_scaler.fit_transform(train_in.iloc[:,1:-1])
scaled_train_out = output_scaler.fit_transform(train_out.iloc[:,1:-1])
scaled_val_in = input_scaler.transform(val_in.iloc[:,1:-1])
scaled_val_out = output_scaler.transform(val_out.iloc[:,1:-1])
scaled_train_in.shape, scaled_train_out.shape, scaled_val_in.shape, scaled_val_out.shape
# Standardizing
input_std_scaler = StandardScaler()
output_std_scaler = StandardScaler()

scaled_std_train_in = input_std_scaler.fit_transform(train_in.iloc[:,1:-1])
scaled_std_train_out = output_std_scaler.fit_transform(train_out.iloc[:,1:-1])
scaled_std_val_in = input_std_scaler.transform(val_in.iloc[:,1:-1])
scaled_std_val_out = output_std_scaler.transform(val_out.iloc[:,1:-1])
scaled_std_train_in.shape, scaled_std_train_out.shape, scaled_std_val_in.shape, scaled_std_val_out.shape<jupyter_output><empty_output><jupyter_text>Baseline model and previous ML model<jupyter_code>%%time
tf.random.set_seed(seed)
mod = dense_neural_net()
history = mod.fit(scaled_train_in, scaled_train_out, validation_data=(scaled_val_in, scaled_val_out), 
                  batch_size=256, epochs=5, verbose=1)
mod.summary()
# get predictions and truth

# pred_nn = pd.DataFrame(output_scaler.inverse_transform(mod.predict(scaled_val_in))).iloc[1:,:]
# true_val = val_out.iloc[:-1,1:4]

pred_nn = pd.DataFrame(output_scaler.inverse_transform(mod.predict(scaled_val_in)))
true_val = val_out.iloc[:,1:-1]

print('Metrics for base model:')
evaluate_mod(true_val, pred_nn)
output_columns = ['Precursor [ug/m3]', 'Gas [ug/m3]', 'Aerosol [ug_m3]']
pred_nn.columns = output_columns

pred_train_nn = pd.DataFrame(output_scaler.inverse_transform(mod.predict(scaled_train_in)))
pred_train_nn.columns = output_columns<jupyter_output>/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:5: UserWarning: Pandas doesn't allow columns to be created via a new attribute name - see https://pandas.pydata.org/pandas-docs/stable/indexing.html#attribute-access
  """<jupyter_text>Baseline model with standardized input.<jupyter_code>%%time
tf.random.set_seed(seed)
mod_std = dense_neural_net()
history = mod_std.fit(scaled_train_in, scaled_train_out, validation_data=(scaled_std_val_in, scaled_std_val_out), 
                  batch_size=256, epochs=5, verbose=1)
mod_std.summary()
# get predictions and truth

# pred_nn = pd.DataFrame(output_std_scaler.inverse_transform(mod.predict(scaled_std_val_in))).iloc[1:,:]
# true_val = val_out.iloc[:-1,1:4]

pred_nn_std = pd.DataFrame(output_std_scaler.inverse_transform(mod_std.predict(scaled_std_val_in)))
true_val = val_out.iloc[:,1:-1]

print('Metrics for base model:')
evaluate_mod(true_val, pred_nn_std)<jupyter_output>Metrics for base model:
RMSE: Precursor: 0.00246, Gas: 0.00952, Aerosols: 0.00884
R2: Precursor: 0.98158, Gas: 0.95517, Aerosols: 0.95286
Hellenger Distance: Precursor: 0.26690, Gas: 0.18135, Aerosols: 0.09445<jupyter_text>Evaluate two interpretation methods for your machine learning solution  Models to evaluate The RNN Model to evaluate:<jupyter_code>from tensorflow.keras.layers import LSTM, SimpleRNN
from tensorflow.keras.models import Sequential

tf.random.set_seed(seed)

# Transform the model for use in RNN
# previously applied transformer: output_std_scaler
rnn_std_train_in = scaled_std_train_in.reshape(-1,1,9)
rnn_std_train_out = scaled_std_train_out.reshape(-1,1,3)
rnn_std_val_in = scaled_std_val_in.reshape(-1,1,9)
rnn_std_val_out = scaled_std_val_out.reshape(-1,1,3)

# Set up the model
mod_rnn_std_2 = Sequential()
mod_rnn_std_2.add(SimpleRNN(128, input_shape=(1, 9), return_sequences=True))
mod_rnn_std_2.add(LSTM(256, return_sequences=True))
mod_rnn_std_2.add(Dense(3))
mod_rnn_std_2.compile(Adam(learning_rate=0.001), "mse")

# fit the model
history = mod_rnn_std_2.fit(rnn_std_train_in, rnn_std_train_out, 
                      validation_data=(rnn_std_val_in, rnn_std_val_out), 
                      batch_size=256, epochs=5, verbose=1, workers=-1)
mod_rnn_std_2.summary()

# get predictions and truth
pred_rnn_std_reshape_2 = mod_rnn_std_2.predict(rnn_std_val_in).reshape(scaled_std_val_out.shape)
pred_rnn_std_2 = pd.DataFrame(output_std_scaler.inverse_transform(pred_rnn_std_reshape_2))
true_val = val_out.iloc[:,1:-1]

print('Metrics for LSTM model:')
evaluate_mod(true_val, pred_rnn_std_2)
output_columns = ['Precursor [ug/m3]', 'Gas [ug/m3]', 'Aerosol [ug_m3]']
pred_rnn_std_2.columns = output_columns

pred_train_rnn_std_reshape_2 = mod_rnn_std_2.predict(rnn_std_train_in).reshape(scaled_std_train_out.shape)
pred_train_rnn_std_2 = pd.DataFrame(output_std_scaler.inverse_transform(pred_train_rnn_std_reshape_2))

pred_train_rnn_std_2.columns = output_columns<jupyter_output><empty_output><jupyter_text>Linear model to evaluate<jupyter_code># from sklearn import linear_model
# lmod = linear_model.LinearRegression().fit(scaled_std_train_in, scaled_std_train_out)

# # get predictions and truth
# pred_lmod_raw = lmod.predict(scaled_std_val_in)
# pred_lmod = pd.DataFrame(output_std_scaler.inverse_transform(pred_rnn_std_reshape_8))
# true_val = val_out.iloc[:,1:-1]

# print('Metrics for Linear model:')
# evaluate_mod(true_val, pred_lmod)
index_columns = ['Time [s]', 'id']
output_columns = ['Precursor [ug/m3]', 'Gas [ug/m3]', 'Aerosol [ug_m3]']
input_columns = list(set(train.columns) - set(index_columns))

def train_model():
    est= {}
    msqe = {}
    for j, col in enumerate(['Precursor [ug/m3]', 'Gas [ug/m3]', 'Aerosol [ug_m3]']):
        e = linear_model.LinearRegression().fit(scaled_std_train_in, scaled_std_train_out[:,j])
#         e = svm.LinearSVR(tol=0.9).fit(train_in_scaled, train_out_scaled[:,j])
        msqe[col] = mean_squared_error(scaled_std_val_out[:,j], e.predict(scaled_std_val_in))
        est[col] = e
    return msqe, est
msqe, est = train_model()

fit_train = {col: est[col].predict(scaled_std_train_in) for col in output_columns}
fit_val = {col: est[col].predict(scaled_std_val_in) for col in output_columns}
train_out_arr = np.concatenate([fit_train[col][:,np.newaxis] for col in output_columns], axis=-1)
val_out_arr = np.concatenate([fit_val[col][:,np.newaxis] for col in output_columns], axis=-1)
linear_mod_train = pd.DataFrame(output_std_scaler.inverse_transform(train_out_arr), columns=output_columns)
linear_mod_val = pd.DataFrame(output_std_scaler.inverse_transform(val_out_arr), columns=output_columns)<jupyter_output><empty_output><jupyter_text>Random Forest Model to evaluate<jupyter_code>from sklearn.ensemble import RandomForestRegressor
# Instantiate model 
model_rf = RandomForestRegressor(n_estimators=20, max_depth=5)
# Train the model with training set
model_rf.fit(scaled_std_train_in, scaled_std_train_out)
# Score
print('Training data score: ' + str(model_rf.score(scaled_std_train_in, scaled_std_train_out)))
print('Testing data score: ' + str(model_rf.score(scaled_std_val_in, scaled_std_val_out)))
# get predictions and truth
pred_rf = pd.DataFrame(output_std_scaler.inverse_transform(model_rf.predict(scaled_std_val_in)))
true_val = val_out.iloc[:,1:4]

print('Metrics for random forest model:')
evaluate_mod(true_val, pred_rf)
output_columns = ['Precursor [ug/m3]', 'Gas [ug/m3]', 'Aerosol [ug_m3]']
pred_rf.columns = output_columns

pred_train_rf = pd.DataFrame(output_std_scaler.inverse_transform(model_rf.predict(scaled_std_train_in)))
pred_train_rf.columns = output_columns<jupyter_output><empty_output><jupyter_text>Simple RNN<jupyter_code>from tensorflow.keras.layers import LSTM, SimpleRNN
from tensorflow.keras.models import Sequential

tf.random.set_seed(seed)

# Transform the model for use in RNN
# previously applied transformer: output_std_scaler
rnn_std_train_in = scaled_std_train_in.reshape(-1,1,9)
rnn_std_train_out = scaled_std_train_out.reshape(-1,1,3)
rnn_std_val_in = scaled_std_val_in.reshape(-1,1,9)
rnn_std_val_out = scaled_std_val_out.reshape(-1,1,3)

# Set up the model
mod_rnn_std_simple = Sequential()
mod_rnn_std_simple.add(SimpleRNN(256, input_shape=(1, 9), return_sequences=True))
mod_rnn_std_simple.add(Dense(3))
mod_rnn_std_simple.compile(Adam(learning_rate=0.001), "mse")

# fit the model
history = mod_rnn_std_simple.fit(rnn_std_train_in, rnn_std_train_out, 
                      validation_data=(rnn_std_val_in, rnn_std_val_out), 
                      batch_size=256, epochs=5, verbose=1, workers=-1)
mod_rnn_std_simple.summary()

# get predictions and truth
pred_rnn_std_reshape_simple = mod_rnn_std_simple.predict(rnn_std_val_in).reshape(scaled_std_val_out.shape)
pred_rnn_std_simple = pd.DataFrame(output_std_scaler.inverse_transform(pred_rnn_std_reshape_simple))
true_val = val_out.iloc[:,1:-1]

print('Metrics for LSTM model:')
evaluate_mod(true_val, pred_rnn_std_simple)
output_columns = ['Precursor [ug/m3]', 'Gas [ug/m3]', 'Aerosol [ug_m3]']
pred_rnn_std_simple.columns = output_columns

pred_train_rnn_std_reshape_simple = mod_rnn_std_simple.predict(rnn_std_train_in).reshape(scaled_std_train_out.shape)
pred_train_rnn_std_simple = pd.DataFrame(output_std_scaler.inverse_transform(pred_train_rnn_std_reshape_simple))

pred_train_rnn_std_simple.columns = output_columns<jupyter_output><empty_output><jupyter_text>LSTM<jupyter_code>from tensorflow.keras.layers import LSTM, SimpleRNN
from tensorflow.keras.models import Sequential

tf.random.set_seed(seed)

# Transform the model for use in RNN
# previously applied transformer: output_std_scaler
rnn_std_train_in = scaled_std_train_in.reshape(-1,1,9)
rnn_std_train_out = scaled_std_train_out.reshape(-1,1,3)
rnn_std_val_in = scaled_std_val_in.reshape(-1,1,9)
rnn_std_val_out = scaled_std_val_out.reshape(-1,1,3)

# Set up the model
mod_rnn_std_lstm = Sequential()
mod_rnn_std_lstm.add(SimpleRNN(256, input_shape=(1, 9), return_sequences=True))
mod_rnn_std_lstm.add(Dense(3))
mod_rnn_std_lstm.compile(Adam(learning_rate=0.001), "mse")

# fit the model
history = mod_rnn_std_lstm.fit(rnn_std_train_in, rnn_std_train_out, 
                      validation_data=(rnn_std_val_in, rnn_std_val_out), 
                      batch_size=256, epochs=5, verbose=1, workers=-1)
mod_rnn_std_lstm.summary()

# get predictions and truth
pred_rnn_std_reshape_lstm = mod_rnn_std_lstm.predict(rnn_std_val_in).reshape(scaled_std_val_out.shape)
pred_rnn_std_lstm = pd.DataFrame(output_std_scaler.inverse_transform(pred_rnn_std_reshape_lstm))
true_val = val_out.iloc[:,1:-1]

print('Metrics for LSTM model:')
evaluate_mod(true_val, pred_rnn_std_lstm)
output_columns = ['Precursor [ug/m3]', 'Gas [ug/m3]', 'Aerosol [ug_m3]']
pred_rnn_std_lstm.columns = output_columns

pred_train_rnn_std_reshape_lstm = mod_rnn_std_lstm.predict(rnn_std_train_in).reshape(scaled_std_train_out.shape)
pred_train_rnn_std_lstm = pd.DataFrame(output_std_scaler.inverse_transform(pred_train_rnn_std_reshape_lstm))

pred_train_rnn_std_lstm.columns = output_columns<jupyter_output><empty_output><jupyter_text>Taylor Diagram<jupyter_code>train_model_outputs = {
    'Linear Train': linear_mod_train,
    'DNN Train': pred_train_nn,
    'RNN Train': pred_train_rnn_std_2,
    'Simple RNN Train': pred_train_rnn_std_simple,
    'LSTM Train': pred_train_rnn_std_lstm,
    'Random Forest Train': pred_train_rf,
    'Linear Val':linear_mod_val,
    'DNN Val': pred_nn,
    'RNN Val': pred_rnn_std_2, 
    'Simple RNN Val': pred_rnn_std_simple,
    'LSTM Val': pred_rnn_std_lstm,
    'Random Forest Val': pred_rf
}
truth = {
    'Train': true_val,
    'Val': true_val
}

def get_stats(truth, models):
    stats = pd.DataFrame(index=['MAE', 'RMSE', 'R', 'Hellinger Distance', 'STDVAR', 'Truth STDVAR', 'Col', 'Model Type'])
    for col in truth.columns.values:
        for model_name, model in models.items():
            mae = (truth[col] - model[col]).abs().mean()
            rmse = np.sqrt(((truth[col] - model[col])**2).mean())
            pearson_corr = ((truth[col]*model[col]).mean() - truth[col].mean()*model[col].mean())/(truth[col].std()*model[col].std())
            hd = hellinger_distance(truth[col].values, model[col].values)
            model_std = model[col].std()
            truth_std = truth[col].std()
            stats[f'{model_name} {col}'] = np.array([mae, rmse, pearson_corr, hd, model_std, truth_std, col, model_name])
    return stats.transpose().reset_index().drop(columns='index') #.reset_index().set_index(['Col', 'Model Type'])

df_all_val = get_stats(truth['Val'], {f'{t} Val': train_model_outputs[f'{t} Val'] for t in ['Linear', 'DNN', 'RNN', 'Simple RNN', 'LSTM', 'Random Forest']})
scores_val = {outname: df_out.set_index('Model Type').drop(columns='Col') for outname, df_out in df_all_val.groupby('Col')}
df_all_train = get_stats(truth['Val'], {f'{t} Train': train_model_outputs[f'{t} Train'] for t in ['Linear', 'DNN', 'RNN', 'Simple RNN', 'LSTM', 'Random Forest']})
scores_train = {outname: df_out.set_index('Model Type').drop(columns='Col') for outname, df_out in df_all_train.groupby('Col')}

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots(figsize=(8,8), subplot_kw=dict(projection='polar'))
corrs = [0.1, 0.25, 0.5, 0.7, 0.9,0.95, 0.99]
theta_corrs = np.arccos(corrs)*180./np.pi
ax.set_thetamax(90)
ax.set_thetagrids(theta_corrs, corrs)
var_space, theta_space = np.meshgrid(np.linspace(0,2,30), np.linspace(0,np.pi/2, 30))
bias_contours = np.sqrt(var_space**2+1-2*var_space*np.cos(theta_space))
contours = ax.contour(theta_space, var_space, bias_contours, colors='k')
ax.clabel(contours, contours.levels, inline=True, fmt='Bias: %2.1f', fontsize=10)
colors = ['r', 'b', 'g', 'c', 'm', 'k']
for modeltype, color in zip(['Linear', 'DNN', 'RNN', 'Simple RNN', 'LSTM', 'Random Forest'], colors):
    r = []
    stdvar_norms = []
    labels = []
    for datatype, dataset in zip(['Train', 'Val'], [scores_train, scores_val]):
        for component in output_columns:
            element_r = float(dataset[component].loc[f'{modeltype} {datatype}', 'R'])
            r.append(element_r)
            std = float(dataset[component].loc[f'{modeltype} {datatype}', 'STDVAR'])
            truth_std = float(dataset[component].loc[f'{modeltype} {datatype}', 'Truth STDVAR'])
            stdvar_norms.append(std/truth_std)
#             ax.text(std/truth_std, np.arccos(element_r),f'{datatype}')
    ax.scatter( stdvar_norms,np.arccos(r), color=color, label=modeltype, s = 25)
ax.legend(loc=2, bbox_to_anchor=(1.05,1))<jupyter_output><empty_output><jupyter_text>Validation Curves Validation curve here was plotted with a different model: this RNN model has its input formatted as 3-d array to show time series.<jupyter_code># Data reprocessing
def time_series(df, colnum):
    '''
    turn dataframe into multi-index dataframe, then to xarray dataset,
    then to xarray dataarray. After that, transpose to the required 
    input shape ([batch, timesteps, feature]).

    returns the 3-d np array above in 2-d (-1,colnum) format so it can be
    standardized or scaled.
    '''
    data_array = df.set_index(['id','Time [s]']).to_xarray().to_array()
    array = data_array.transpose('id','Time [s]', 'variable').values

    return array.reshape(-1, colnum) # reshaped to scaled with xx_std_scaler

# np.shape(time_series(train_in, 9))
# # (2014600, 9)

# actual shape is (1400, 1439, 9)
# Data preprocessing
# Standardizing
input_std_scaler = StandardScaler()
output_std_scaler = StandardScaler()

scaled_std_train_in_ts = input_std_scaler.fit_transform(time_series(train_in, 9))
scaled_std_train_out_ts = output_std_scaler.fit_transform(time_series(train_out, 3))
scaled_std_val_in_ts = input_std_scaler.transform(time_series(val_in, 9))
scaled_std_val_out_ts = output_std_scaler.transform(time_series(val_out, 3))

# reshape. note this time data is in [batch, timesteps, feature] size.
rnn_std_train_in_ts = scaled_std_train_in_ts.reshape(-1,1439,9)
rnn_std_train_out_ts = scaled_std_train_out_ts.reshape(-1,1439,3)
rnn_std_val_in_ts = scaled_std_val_in_ts.reshape(-1,1439,9)
rnn_std_val_out_ts = scaled_std_val_out_ts.reshape(-1,1439,3)
from tensorflow.keras.layers import LSTM, SimpleRNN
from tensorflow.keras.models import Sequential

tf.random.set_seed(seed)

# Set up the model
mod_rnn_std_8 = Sequential()
mod_rnn_std_8.add(SimpleRNN(256, input_shape=(1439, 9), return_sequences=True))
# mod_rnn_std_8.add(LSTM(256, return_sequences=True))
mod_rnn_std_8.add(Dense(3))
mod_rnn_std_8.compile(Adam(learning_rate=0.001), "mse")

# fit the model
# history = mod_rnn_std_8.fit(rnn_std_train_in_ts, rnn_std_train_out_ts, 
#                       validation_data=(rnn_std_val_in_ts, rnn_std_val_out_ts), 
#                       batch_size=10, epochs=5, verbose=1, workers=-1)

history = mod_rnn_std_8.fit(rnn_std_train_in_ts, rnn_std_train_out_ts, 
                      validation_data=(rnn_std_val_in_ts, rnn_std_val_out_ts), 
                      batch_size=256, epochs=1000, verbose=1, workers=-1)

mod_rnn_std_8.summary()

# get predictions and truth
pred_rnn_std_8 = mod_rnn_std_8.predict(rnn_std_val_in_ts).reshape(-1,3)
pred_rnn_std_reshape_8 = mod_rnn_std_8.predict(rnn_std_val_in_ts).reshape(scaled_std_val_out.shape)
pred_rnn_std_reverse_8 = pd.DataFrame(output_std_scaler.inverse_transform(pred_rnn_std_reshape_8)).iloc[1:,:]
true_val = val_out.iloc[:-1,1:4]

print('Metrics for LSTM model:')
evaluate_mod(true_val, pred_rnn_std_reverse_8)
def plot_history():
    plt.figure()
    plt.xlabel('Epoch')
    plt.ylabel('Mean Square Error')
    plt.plot(hist['epoch'], hist['loss'], label='Train Error')
    plt.plot(hist['epoch'], hist['val_loss'], label = 'Val Error')
    plt.legend()

hist = pd.DataFrame(history.history)
hist['epoch'] = history.epoch

plot_history()
def plot_history_short():
    plt.figure()
    plt.xlabel('Epoch')
    plt.ylabel('Mean Square Error')
    plt.plot(hist['epoch'], hist['loss'], label='Train Error')
    plt.plot(hist['epoch'], hist['val_loss'], label = 'Val Error')
    plt.legend()
    plt.xlim([0,100])
    plt.ylim([0,0.10])

plot_history_short()<jupyter_output><empty_output><jupyter_text>Compare interpretation of baseline with your approach<jupyter_code>j=0
print('====================================')
print(output_columns[j])
print('====================================')
foo = scores_val[output_columns[j]]
for col in foo.columns:
    foo[col] = pd.to_numeric(foo[col])
foo.style.background_gradient(cmap='Reds', subset=['MAE', 'RMSE', 'Hellinger Distance']).background_gradient(cmap='Greens', subset=['R'])
j=0
print('====================================')
print(output_columns[j])
print('====================================')
foo = scores_train[output_columns[j]]
for col in foo.columns:
    foo[col] = pd.to_numeric(foo[col])
foo.style.background_gradient(cmap='Reds', subset=['MAE', 'RMSE', 'Hellinger Distance']).background_gradient(cmap='Greens', subset=['R'])<jupyter_output>====================================
Precursor [ug/m3]
====================================
{'fraction_non_alphanumeric': 0.13671432980146075, 'fraction_numerical': 0.017830813016493502, 'mean_word_length': 5.112764619576609, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 51, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '254282', 'license': 'MIT', 'n_tokens_mistral': 10982, 'n_tokens_neox': 10607, 'n_words': 1952, 'path': 'gecko_team_17_thursday.ipynb', 'repo_name': 'ShutingZhai/ai4ess-hackathon-2020', 'size': '29163'}