
%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 15:54:32 2023


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


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


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


data = df_bottle_ta[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_ta[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

((4151, 4), (3113, 4), (1038, 4), (3113, 1), (1038, 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 15:54:55.614640: 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.

78/78 - 1s - loss: 5058152.5000 - val_loss: 4826869.5000 - 1s/epoch - 16ms/step
Epoch 2/700
78/78 - 0s - loss: 5024779.0000 - val_loss: 4675571.0000 - 209ms/epoch - 3ms/step
Epoch 3/700
78/78 - 0s - loss: 4959481.0000 - val_loss: 4727107.0000 - 220ms/epoch - 3ms/step
Epoch 4/700
78/78 - 0s - loss: 4857733.5000 - val_loss: 4667107.5000 - 214ms/epoch - 3ms/step
Epoch 5/700
78/78 - 0s - loss: 4721133.0000 - val_loss: 4553366.0000 - 216ms/epoch - 3ms/step
Epoch 6/700
78/78 - 0s - loss: 4552999.0000 - val_loss: 4386781.0000 - 214ms/epoch - 3ms/step
Epoch 7/700
78/78 - 0s - loss: 4357344.0000 - val_loss: 4223865.0000 - 223ms/epoch - 3ms/step
Epoch 8/700
78/78 - 0s - loss: 4138198.7500 - val_loss: 4003955.5000 - 223ms/epoch - 3ms/step
Epoch 9/700
78/78 - 0s - loss: 3899918.0000 - val_loss: 3768097.2500 - 220ms/epoch - 3ms/step
Epoch 10/700
78/78 - 0s - loss: 3646812.7500 - val_loss: 3509264.7500 - 205ms/epoch - 3ms/step
Epoch 11/700
78/78 - 0s - loss: 3383293.7500 - val_loss: 3255060.7500 - 201ms/epoch - 3ms/step
Epoch 12/700
78/78 - 0s - loss: 3113344.5000 - val_loss: 2984986.2500 - 212ms/epoch - 3ms/step
Epoch 13/700
78/78 - 0s - loss: 2841055.0000 - val_loss: 2700406.5000 - 213ms/epoch - 3ms/step
Epoch 14/700
78/78 - 0s - loss: 2570210.2500 - val_loss: 2430672.7500 - 202ms/epoch - 3ms/step
Epoch 15/700
78/78 - 0s - loss: 2304376.0000 - val_loss: 2175756.7500 - 196ms/epoch - 3ms/step
Epoch 16/700
78/78 - 0s - loss: 2046803.5000 - val_loss: 1916426.1250 - 196ms/epoch - 3ms/step
Epoch 17/700
78/78 - 0s - loss: 1800169.1250 - val_loss: 1680347.3750 - 207ms/epoch - 3ms/step
Epoch 18/700
78/78 - 0s - loss: 1566914.8750 - val_loss: 1452284.3750 - 207ms/epoch - 3ms/step
Epoch 19/700
78/78 - 0s - loss: 1349124.3750 - val_loss: 1244729.6250 - 206ms/epoch - 3ms/step
Epoch 20/700
78/78 - 0s - loss: 1148267.7500 - val_loss: 1050262.6250 - 207ms/epoch - 3ms/step
Epoch 21/700
78/78 - 0s - loss: 965429.0625 - val_loss: 877454.5625 - 204ms/epoch - 3ms/step
Epoch 22/700
78/78 - 0s - loss: 801303.1250 - val_loss: 725895.9375 - 210ms/epoch - 3ms/step
Epoch 23/700
78/78 - 0s - loss: 656018.6875 - val_loss: 590227.0000 - 207ms/epoch - 3ms/step
Epoch 24/700
78/78 - 0s - loss: 529350.8750 - val_loss: 473904.5000 - 209ms/epoch - 3ms/step
Epoch 25/700
78/78 - 0s - loss: 420552.9062 - val_loss: 370960.0000 - 210ms/epoch - 3ms/step
Epoch 26/700
78/78 - 0s - loss: 328881.5938 - val_loss: 289870.0000 - 209ms/epoch - 3ms/step
Epoch 27/700
78/78 - 0s - loss: 252752.5781 - val_loss: 219759.5625 - 208ms/epoch - 3ms/step
Epoch 28/700
78/78 - 0s - loss: 190713.0000 - val_loss: 162043.5312 - 204ms/epoch - 3ms/step
Epoch 29/700
78/78 - 0s - loss: 141315.9219 - val_loss: 118471.8594 - 212ms/epoch - 3ms/step
Epoch 30/700
78/78 - 0s - loss: 102638.0625 - val_loss: 85916.8672 - 202ms/epoch - 3ms/step
Epoch 31/700
78/78 - 0s - loss: 73053.1562 - val_loss: 59968.0391 - 197ms/epoch - 3ms/step
Epoch 32/700
78/78 - 0s - loss: 50931.4570 - val_loss: 40789.4688 - 203ms/epoch - 3ms/step
Epoch 33/700
78/78 - 0s - loss: 34799.3828 - val_loss: 26687.8867 - 200ms/epoch - 3ms/step
Epoch 34/700
78/78 - 0s - loss: 23261.2969 - val_loss: 17565.3340 - 206ms/epoch - 3ms/step
Epoch 35/700
78/78 - 0s - loss: 15288.5400 - val_loss: 11896.5439 - 206ms/epoch - 3ms/step
Epoch 36/700
78/78 - 0s - loss: 9826.0518 - val_loss: 7006.7466 - 208ms/epoch - 3ms/step
Epoch 37/700
78/78 - 0s - loss: 6190.0757 - val_loss: 3988.5349 - 209ms/epoch - 3ms/step
Epoch 38/700
78/78 - 0s - loss: 3915.7104 - val_loss: 2430.5886 - 213ms/epoch - 3ms/step
Epoch 39/700
78/78 - 0s - loss: 2438.5886 - val_loss: 1863.3517 - 216ms/epoch - 3ms/step
Epoch 40/700
78/78 - 0s - loss: 1586.0765 - val_loss: 1052.0110 - 209ms/epoch - 3ms/step
Epoch 41/700
78/78 - 0s - loss: 1066.4336 - val_loss: 667.0111 - 206ms/epoch - 3ms/step
Epoch 42/700
78/78 - 0s - loss: 695.5529 - val_loss: 415.7258 - 207ms/epoch - 3ms/step
Epoch 43/700
78/78 - 0s - loss: 541.9031 - val_loss: 240.6377 - 207ms/epoch - 3ms/step
Epoch 44/700
78/78 - 0s - loss: 421.1344 - val_loss: 244.1670 - 201ms/epoch - 3ms/step
Epoch 45/700
78/78 - 0s - loss: 370.2583 - val_loss: 145.9110 - 207ms/epoch - 3ms/step
Epoch 46/700
78/78 - 0s - loss: 352.1632 - val_loss: 135.6085 - 205ms/epoch - 3ms/step
Epoch 47/700
78/78 - 0s - loss: 393.1107 - val_loss: 148.2265 - 200ms/epoch - 3ms/step
Epoch 48/700
78/78 - 0s - loss: 341.7621 - val_loss: 130.0451 - 213ms/epoch - 3ms/step
Epoch 49/700
78/78 - 0s - loss: 313.5030 - val_loss: 138.7267 - 208ms/epoch - 3ms/step
Epoch 50/700
78/78 - 0s - loss: 388.6017 - val_loss: 132.5679 - 201ms/epoch - 3ms/step
Epoch 51/700
78/78 - 0s - loss: 408.2587 - val_loss: 124.7427 - 201ms/epoch - 3ms/step
Epoch 52/700
78/78 - 0s - loss: 348.2925 - val_loss: 124.7329 - 209ms/epoch - 3ms/step
Epoch 53/700
78/78 - 0s - loss: 386.9626 - val_loss: 136.0415 - 210ms/epoch - 3ms/step
Epoch 54/700
78/78 - 0s - loss: 315.8572 - val_loss: 180.2669 - 202ms/epoch - 3ms/step
Epoch 55/700
78/78 - 0s - loss: 343.1501 - val_loss: 145.8454 - 201ms/epoch - 3ms/step
Epoch 56/700
78/78 - 0s - loss: 382.6131 - val_loss: 169.8068 - 201ms/epoch - 3ms/step
Epoch 57/700
78/78 - 0s - loss: 337.7169 - val_loss: 131.2766 - 200ms/epoch - 3ms/step
Epoch 58/700
78/78 - 0s - loss: 298.1558 - val_loss: 124.1429 - 211ms/epoch - 3ms/step
Epoch 59/700
78/78 - 0s - loss: 414.9820 - val_loss: 145.4434 - 206ms/epoch - 3ms/step
Epoch 60/700
78/78 - 0s - loss: 366.7910 - val_loss: 138.3793 - 201ms/epoch - 3ms/step
Epoch 61/700
78/78 - 0s - loss: 326.3562 - val_loss: 207.3313 - 201ms/epoch - 3ms/step
Epoch 62/700
78/78 - 0s - loss: 393.8535 - val_loss: 131.2289 - 199ms/epoch - 3ms/step
Epoch 63/700
78/78 - 0s - loss: 372.0420 - val_loss: 142.0009 - 204ms/epoch - 3ms/step
Epoch 64/700
78/78 - 0s - loss: 355.8058 - val_loss: 131.4643 - 205ms/epoch - 3ms/step
Epoch 65/700
78/78 - 0s - loss: 303.9357 - val_loss: 144.1300 - 210ms/epoch - 3ms/step
Epoch 66/700
78/78 - 0s - loss: 374.0136 - val_loss: 128.3093 - 213ms/epoch - 3ms/step
Epoch 67/700
78/78 - 0s - loss: 372.1476 - val_loss: 157.0159 - 209ms/epoch - 3ms/step
Epoch 68/700
78/78 - 0s - loss: 389.1917 - val_loss: 141.8057 - 215ms/epoch - 3ms/step
Epoch 69/700
78/78 - 0s - loss: 334.5516 - val_loss: 124.6738 - 212ms/epoch - 3ms/step
Epoch 70/700
78/78 - 0s - loss: 309.6681 - val_loss: 115.7426 - 209ms/epoch - 3ms/step
Epoch 71/700
78/78 - 0s - loss: 335.8249 - val_loss: 116.5471 - 202ms/epoch - 3ms/step
Epoch 72/700
78/78 - 0s - loss: 306.6313 - val_loss: 144.6857 - 210ms/epoch - 3ms/step
Epoch 73/700
78/78 - 0s - loss: 341.2404 - val_loss: 149.5908 - 210ms/epoch - 3ms/step
Epoch 74/700
78/78 - 0s - loss: 347.8168 - val_loss: 118.0576 - 205ms/epoch - 3ms/step
Epoch 75/700
78/78 - 0s - loss: 377.9933 - val_loss: 128.6920 - 204ms/epoch - 3ms/step
Epoch 76/700
78/78 - 0s - loss: 295.1218 - val_loss: 137.8609 - 204ms/epoch - 3ms/step
Epoch 77/700
78/78 - 0s - loss: 315.6928 - val_loss: 126.2048 - 205ms/epoch - 3ms/step
Epoch 78/700
78/78 - 0s - loss: 378.3812 - val_loss: 125.3822 - 206ms/epoch - 3ms/step
Epoch 79/700
78/78 - 0s - loss: 305.4313 - val_loss: 169.6523 - 201ms/epoch - 3ms/step
Epoch 80/700
78/78 - 0s - loss: 380.0275 - val_loss: 164.7089 - 205ms/epoch - 3ms/step
Epoch 81/700
78/78 - 0s - loss: 367.9416 - val_loss: 120.7348 - 204ms/epoch - 3ms/step
Epoch 82/700
78/78 - 0s - loss: 295.1768 - val_loss: 150.8757 - 200ms/epoch - 3ms/step
Epoch 83/700
78/78 - 0s - loss: 336.2103 - val_loss: 135.8440 - 206ms/epoch - 3ms/step
Epoch 84/700
78/78 - 0s - loss: 295.1935 - val_loss: 119.2881 - 205ms/epoch - 3ms/step
Epoch 85/700
78/78 - 0s - loss: 355.3746 - val_loss: 119.9988 - 203ms/epoch - 3ms/step
Epoch 86/700
78/78 - 0s - loss: 393.3821 - val_loss: 149.6162 - 207ms/epoch - 3ms/step
Epoch 87/700
78/78 - 0s - loss: 356.8006 - val_loss: 150.9677 - 210ms/epoch - 3ms/step
Epoch 88/700
78/78 - 0s - loss: 343.5536 - val_loss: 122.6854 - 211ms/epoch - 3ms/step
Epoch 89/700
78/78 - 0s - loss: 289.3266 - val_loss: 117.6777 - 201ms/epoch - 3ms/step
Epoch 90/700
78/78 - 0s - loss: 277.4659 - val_loss: 127.4799 - 215ms/epoch - 3ms/step


model.save('models/nn_regression_ta_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_ta_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 = 104.50
MSE on test set     = 107.46

R squared on training set = 0.982
R squared on test set     = 0.983


fig, ax = plt.subplots(figsize=(6,6))
_ = sns.scatterplot(x=y_test.ravel(), y=y_pred_test.ravel(), ax=ax)
_ = ax.set(xlabel='observed TA', ylabel='predicted TA', 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 + ['TA observed', 'TA predicted'])
df_test['TA residuals'] = df_test['TA observed'] - df_test['TA predicted']
df_test.to_csv('results/bottle_data_test_ta_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.973
Fold 2 test set R squared: 0.983
Fold 3 test set R squared: 0.982
Fold 4 test set R squared: 0.981
Fold 5 test set R squared: 0.980

Best R squared:  0.983
Worst R squared: 0.973
Mean R squared:  0.980

Neural network model to predict TA¶

Load TA 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¶