
%matplotlib inline
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import os, datetime, warnings
print('Last updated on {}'.format(datetime.datetime.now().ctime()))

Last updated on Sun Feb  5 13:24:53 2023


import sns_settings
sns.set_context('paper')
pd.options.display.max_columns = 50
warnings.filterwarnings('ignore')


df_bottle_dic = pd.read_csv('data/bottle_data_DIC_prepared.csv', parse_dates=['Date'],
                            index_col=0, na_values=['<undefined>',-9999.])
df_bottle_dic['log_Chl'] = np.log(df_bottle_dic.Chl)
df_bottle_dic['log_KD490'] = np.log(df_bottle_dic.KD490)


features = ['Depth', 'Temperature', 'Salinity', 'pCO2_atm']
target = ['DIC']
varlist = features + target
fg = sns.pairplot(df_bottle_dic, vars=varlist, hue='Season', plot_kws={'alpha':0.7}, diag_kind='hist')


data = df_bottle_dic[varlist]
corr_mat = data.corr()
fig, ax = plt.subplots(figsize=(7,7))
_ = sns.heatmap(corr_mat, ax=ax, cmap='vlag', center=0, square=True, annot=True, annot_kws={'fontsize':9})
_ = ax.set_title('Correlation')


from sklearn.model_selection import train_test_split, cross_val_score

data = df_bottle_dic[features + target + ['Season']].dropna()

X = data[features].values
y = data[target].values

X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=data.Season.values, random_state=77)
X.shape, X_train.shape, X_test.shape, y_train.shape, y_test.shape

((4350, 4), (3262, 4), (1088, 4), (3262, 1), (1088, 1))


import tensorflow as tf
from tensorflow import keras

keras.utils.set_random_seed(42) # make things reproducible
n_hidden = 256 # number of nodes in hidden layers
alpha=0.01

model = keras.models.Sequential([
    keras.layers.BatchNormalization(),
    keras.layers.Dense(n_hidden, input_shape=X_train.shape[1:]),
    keras.layers.LeakyReLU(alpha=alpha),
    keras.layers.BatchNormalization(),
    keras.layers.Dense(n_hidden),
    keras.layers.LeakyReLU(alpha=alpha),
    keras.layers.BatchNormalization(),
    keras.layers.Dense(y_train.shape[1])
])

early_stopping_cb = keras.callbacks.EarlyStopping(patience=20, restore_best_weights=True)
model.compile(loss='mean_squared_error', optimizer=keras.optimizers.Adam())
history = model.fit(X_train, y_train, epochs=700, verbose=2, validation_split=0.2, callbacks=[early_stopping_cb])

Epoch 1/700

2023-02-05 13:25:16.495891: I tensorflow/core/platform/cpu_feature_guard.cc:151] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  SSE4.1 SSE4.2 AVX AVX2 FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.

82/82 - 1s - loss: 4295505.5000 - val_loss: 4163096.7500 - 1s/epoch - 15ms/step
Epoch 2/700
82/82 - 0s - loss: 4261746.5000 - val_loss: 4100688.7500 - 222ms/epoch - 3ms/step
Epoch 3/700
82/82 - 0s - loss: 4195334.0000 - val_loss: 4025331.0000 - 228ms/epoch - 3ms/step
Epoch 4/700
82/82 - 0s - loss: 4091928.7500 - val_loss: 3923033.7500 - 228ms/epoch - 3ms/step
Epoch 5/700
82/82 - 0s - loss: 3953896.2500 - val_loss: 3800570.5000 - 234ms/epoch - 3ms/step
Epoch 6/700
82/82 - 0s - loss: 3785024.0000 - val_loss: 3641637.7500 - 221ms/epoch - 3ms/step
Epoch 7/700
82/82 - 0s - loss: 3590033.0000 - val_loss: 3446261.7500 - 212ms/epoch - 3ms/step
Epoch 8/700
82/82 - 0s - loss: 3373499.5000 - val_loss: 3227750.5000 - 215ms/epoch - 3ms/step
Epoch 9/700
82/82 - 0s - loss: 3140517.2500 - val_loss: 2988305.7500 - 214ms/epoch - 3ms/step
Epoch 10/700
82/82 - 0s - loss: 2895695.0000 - val_loss: 2754291.0000 - 209ms/epoch - 3ms/step
Epoch 11/700
82/82 - 0s - loss: 2643959.7500 - val_loss: 2503563.7500 - 215ms/epoch - 3ms/step
Epoch 12/700
82/82 - 0s - loss: 2389690.7500 - val_loss: 2241687.0000 - 217ms/epoch - 3ms/step
Epoch 13/700
82/82 - 0s - loss: 2137393.5000 - val_loss: 1994962.5000 - 228ms/epoch - 3ms/step
Epoch 14/700
82/82 - 0s - loss: 1890703.3750 - val_loss: 1747408.3750 - 219ms/epoch - 3ms/step
Epoch 15/700
82/82 - 0s - loss: 1653496.2500 - val_loss: 1512308.5000 - 206ms/epoch - 3ms/step
Epoch 16/700
82/82 - 0s - loss: 1428506.7500 - val_loss: 1297011.5000 - 210ms/epoch - 3ms/step
Epoch 17/700
82/82 - 0s - loss: 1218496.0000 - val_loss: 1097465.7500 - 214ms/epoch - 3ms/step
Epoch 18/700
82/82 - 0s - loss: 1025240.2500 - val_loss: 923995.0625 - 210ms/epoch - 3ms/step
Epoch 19/700
82/82 - 0s - loss: 850363.0000 - val_loss: 752726.3750 - 253ms/epoch - 3ms/step
Epoch 20/700
82/82 - 0s - loss: 694582.5000 - val_loss: 605740.5625 - 226ms/epoch - 3ms/step
Epoch 21/700
82/82 - 0s - loss: 558044.8750 - val_loss: 481607.5312 - 218ms/epoch - 3ms/step
Epoch 22/700
82/82 - 0s - loss: 440643.9062 - val_loss: 375302.5312 - 219ms/epoch - 3ms/step
Epoch 23/700
82/82 - 0s - loss: 341732.6875 - val_loss: 287770.6875 - 209ms/epoch - 3ms/step
Epoch 24/700
82/82 - 0s - loss: 259785.8594 - val_loss: 215465.5781 - 208ms/epoch - 3ms/step
Epoch 25/700
82/82 - 0s - loss: 193552.6250 - val_loss: 159422.8594 - 212ms/epoch - 3ms/step
Epoch 26/700
82/82 - 0s - loss: 141151.8750 - val_loss: 113300.4609 - 211ms/epoch - 3ms/step
Epoch 27/700
82/82 - 0s - loss: 100695.7266 - val_loss: 78104.9531 - 211ms/epoch - 3ms/step
Epoch 28/700
82/82 - 0s - loss: 70128.2344 - val_loss: 52440.6914 - 215ms/epoch - 3ms/step
Epoch 29/700
82/82 - 0s - loss: 47845.6094 - val_loss: 35710.0391 - 218ms/epoch - 3ms/step
Epoch 30/700
82/82 - 0s - loss: 31863.4727 - val_loss: 22918.5586 - 217ms/epoch - 3ms/step
Epoch 31/700
82/82 - 0s - loss: 20791.2988 - val_loss: 14663.5146 - 228ms/epoch - 3ms/step
Epoch 32/700
82/82 - 0s - loss: 13323.1104 - val_loss: 8602.0479 - 224ms/epoch - 3ms/step
Epoch 33/700
82/82 - 0s - loss: 8411.1152 - val_loss: 5262.0317 - 223ms/epoch - 3ms/step
Epoch 34/700
82/82 - 0s - loss: 5299.3701 - val_loss: 3135.6426 - 220ms/epoch - 3ms/step
Epoch 35/700
82/82 - 0s - loss: 3411.6531 - val_loss: 1845.9076 - 228ms/epoch - 3ms/step
Epoch 36/700
82/82 - 0s - loss: 2224.8977 - val_loss: 1297.6614 - 228ms/epoch - 3ms/step
Epoch 37/700
82/82 - 0s - loss: 1614.6608 - val_loss: 818.2849 - 217ms/epoch - 3ms/step
Epoch 38/700
82/82 - 0s - loss: 1242.4725 - val_loss: 715.1530 - 221ms/epoch - 3ms/step
Epoch 39/700
82/82 - 0s - loss: 978.2239 - val_loss: 774.4937 - 220ms/epoch - 3ms/step
Epoch 40/700
82/82 - 0s - loss: 876.8972 - val_loss: 662.7278 - 226ms/epoch - 3ms/step
Epoch 41/700
82/82 - 0s - loss: 802.7866 - val_loss: 704.6508 - 221ms/epoch - 3ms/step
Epoch 42/700
82/82 - 0s - loss: 804.9638 - val_loss: 745.8519 - 229ms/epoch - 3ms/step
Epoch 43/700
82/82 - 0s - loss: 829.6318 - val_loss: 789.5859 - 222ms/epoch - 3ms/step
Epoch 44/700
82/82 - 0s - loss: 783.4631 - val_loss: 631.9388 - 227ms/epoch - 3ms/step
Epoch 45/700
82/82 - 0s - loss: 741.9945 - val_loss: 679.4446 - 223ms/epoch - 3ms/step
Epoch 46/700
82/82 - 0s - loss: 794.3704 - val_loss: 822.0247 - 216ms/epoch - 3ms/step
Epoch 47/700
82/82 - 0s - loss: 733.0338 - val_loss: 822.8168 - 220ms/epoch - 3ms/step
Epoch 48/700
82/82 - 0s - loss: 776.0649 - val_loss: 904.6737 - 217ms/epoch - 3ms/step
Epoch 49/700
82/82 - 0s - loss: 802.7091 - val_loss: 829.1224 - 222ms/epoch - 3ms/step
Epoch 50/700
82/82 - 0s - loss: 774.8091 - val_loss: 932.9916 - 213ms/epoch - 3ms/step
Epoch 51/700
82/82 - 0s - loss: 841.3947 - val_loss: 715.3143 - 225ms/epoch - 3ms/step
Epoch 52/700
82/82 - 0s - loss: 747.3798 - val_loss: 773.1365 - 225ms/epoch - 3ms/step
Epoch 53/700
82/82 - 0s - loss: 764.3210 - val_loss: 717.6651 - 215ms/epoch - 3ms/step
Epoch 54/700
82/82 - 0s - loss: 766.6048 - val_loss: 861.3317 - 220ms/epoch - 3ms/step
Epoch 55/700
82/82 - 0s - loss: 793.9297 - val_loss: 939.0875 - 215ms/epoch - 3ms/step
Epoch 56/700
82/82 - 0s - loss: 805.8483 - val_loss: 832.5061 - 222ms/epoch - 3ms/step
Epoch 57/700
82/82 - 0s - loss: 745.0978 - val_loss: 764.8958 - 222ms/epoch - 3ms/step
Epoch 58/700
82/82 - 0s - loss: 712.9371 - val_loss: 727.1599 - 217ms/epoch - 3ms/step
Epoch 59/700
82/82 - 0s - loss: 796.3814 - val_loss: 810.2433 - 208ms/epoch - 3ms/step
Epoch 60/700
82/82 - 0s - loss: 717.5599 - val_loss: 881.5594 - 210ms/epoch - 3ms/step
Epoch 61/700
82/82 - 0s - loss: 746.7147 - val_loss: 736.1664 - 209ms/epoch - 3ms/step
Epoch 62/700
82/82 - 0s - loss: 767.3443 - val_loss: 871.4409 - 204ms/epoch - 2ms/step
Epoch 63/700
82/82 - 0s - loss: 756.6715 - val_loss: 785.0773 - 210ms/epoch - 3ms/step
Epoch 64/700
82/82 - 0s - loss: 770.8275 - val_loss: 799.9147 - 213ms/epoch - 3ms/step


model.save('models/nn_regression_dic_nosat.h5')


df_history = pd.DataFrame(history.history)
df_history.index.name = 'epoch'
df_history = df_history.reset_index()
df_history.to_csv('results/nn_regression_history_dic_nosat.csv')

fig, ax = plt.subplots(figsize=(6, 6))
_ = sns.lineplot(x=df_history.epoch-0.5, y='loss', data=df_history, ax=ax, label='training set')
_ = sns.lineplot(x=df_history.epoch, y='val_loss', data=df_history, ax=ax, label='validation set')
_ = ax.set(ylabel = 'MSE')
# _ = ax.set(yscale='log')


from sklearn.metrics import r2_score

print('MSE on training set = {:.2f}'.format(model.evaluate(X_train, y_train, verbose=0)))
print('MSE on test set     = {:.2f}\n'.format(model.evaluate(X_test, y_test, verbose=0)))

y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)

print('R squared on training set = {:.3f}'.format(r2_score(y_train, y_pred_train)))
print('R squared on test set     = {:.3f}'.format(r2_score(y_test, y_pred_test)))

MSE on training set = 560.60
MSE on test set     = 569.15

R squared on training set = 0.913
R squared on test set     = 0.913


fig, ax = plt.subplots(figsize=(6,6))
_ = sns.scatterplot(x=y_test.ravel(), y=y_pred_test.ravel(), ax=ax)
_ = ax.set(xlabel='observed DIC', ylabel='predicted DIC', title='Test dataset')
_ = ax.axis('equal')


# save test set features, target & predictions
df_test = pd.DataFrame(np.c_[X_test, y_test, y_pred_test], columns = features + ['DIC observed', 'DIC predicted'])
df_test['DIC residuals'] = df_test['DIC observed'] - df_test['DIC predicted']
df_test.to_csv('results/bottle_data_test_dic_nosat.csv')


from sklearn.model_selection import KFold

kf = KFold(n_splits=5, shuffle=True, random_state=42)
score_vals = [] # store score values

nn_reg = keras.models.Sequential([
    keras.layers.BatchNormalization(),
    keras.layers.Dense(n_hidden, input_shape=X_train.shape[1:]),
    keras.layers.LeakyReLU(alpha=alpha),
    keras.layers.BatchNormalization(),
    keras.layers.Dense(n_hidden),
    keras.layers.LeakyReLU(alpha=alpha),
    keras.layers.BatchNormalization(),
    keras.layers.Dense(y_train.shape[1])
])
nn_reg.compile(loss='mean_squared_error', optimizer=keras.optimizers.Adam())

for k, (train_idx, test_idx) in enumerate(kf.split(X_train)):
    X_tr, X_te = X_train[train_idx], X_train[test_idx]
    y_tr, y_te = y_train[train_idx], y_train[test_idx]
    history_cv = nn_reg.fit(X_tr, y_tr, epochs=700, verbose=0, validation_split=0.2, callbacks=[early_stopping_cb])
    y_pred = nn_reg.predict(X_te)
    score = r2_score(y_te, y_pred)
    score_vals.append(score)
    print('Fold {} test set R squared: {:.3f}'.format(k+1, score))

scores = np.array(score_vals)
print('\nBest R squared:  {:.3f}'.format(scores.max()))
print('Worst R squared: {:.3f}'.format(scores.min()))
print('Mean R squared:  {:.3f}'.format(scores.mean()))

Fold 1 test set R squared: 0.909
Fold 2 test set R squared: 0.914
Fold 3 test set R squared: 0.905
Fold 4 test set R squared: 0.889
Fold 5 test set R squared: 0.903

Best R squared:  0.914
Worst R squared: 0.889
Mean R squared:  0.904

Neural network model to predict DIC¶

Load DIC bottle data¶

Select input features and target variable¶

Split data into training and test sets¶

Train Neural Network regression¶

Save trained model¶

Learning curve¶

MSE & $R^2$¶

K-fold cross-validation¶