Bilateral LSTM e GRU con dati di estensimetro

Avevo gia' provato ad usare LSTM per dati derivanti da un estensimetro qui ... adesso ci riprovo con la variante Bilateral LSTM e GRU prendendo spunto da questa pagina 

Questi sono i dati raw. I dati sono campionati ogni 20 minuti e sono presenti diverse fasi di un movimento di versante

Prima di iniziare un paio di considerazioni

In questo grafico sono riassunti gli spostamenti registrati dall'estensimetro, la velocita' di movimento ed i dati pluviometrici

Dalla cross correlazione tra il dato meteo e quello di spostamento si un lag di 27 giorni

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

df_b = pd.read_csv('est2.csv', sep=',',header=0, parse_dates = ['Date'], date_format = '%Y-%m-%d %H:%M:%S')

pluvio = pd.read_csv('arpa.csv', sep=';',header=0, parse_dates = ['Date'], date_format = '%d/%m/%Y %H:%M:%S')

df_b_res = df_b.resample('24h', on='Date').mean()
pluvio_res = pluvio.resample('24h', on='Date').mean()

print(df_b_res)
print(pluvio)

df_b_res['Diff'] = df_b_res['Est'].diff()
df_b_res['Inv'] = 1/df_b_res['Diff']

time = pd.date_range('2024-02-01', periods=214, freq='D')
print(time)
print(pluvio_res.iloc[:,0])
print(df_b_res.iloc[:,0])

df = pd.DataFrame({'series1': df_b_res.iloc[:,0], 'series2': pluvio_res.iloc[:,0]}, index=time)
cross_corr = df['series1'].corr(df['series2'])
def cross_correlation(series1, series2, max_lag):
    corr = []
    for lag in range(-max_lag, max_lag + 1):
        corr.append(df['series1'].shift(lag).corr(df['series2']))
    return corr
max_lag = 40  # You can adjust the maximum lag
lags = range(-max_lag, max_lag + 1)
correlation_values = cross_correlation(df['series1'], df['series2'], max_lag)

# Find the lag that maximizes the correlation
best_lag = lags[np.argmax(correlation_values)]
print(f"The best lag is: {best_lag} days")

#print(df_b_res)

plt.figure(figsize=(11,7))
plt.plot(df_b_res['Est'],"-b",label="Est. (0.1 mm)")
plt.plot(df_b_res['Diff']*35,"-r",label="35*Vel. (0.1 mm/day)")
plt.plot(pluvio_res['Prec'],"-g",label="Rain (mm)")


plt.title('Velocita')
plt.xlabel('Date')
plt.ylabel('Estensimetro 0.1mm')
plt.legend(loc="lower left")

plt.show()

plt.figure(figsize=(11,7))
plt.plot(df_b_res['Inv'],"+b",label="Est.")
plt.title('1/Velocita')
plt.xlabel('Date')
plt.ylabel('1/V day/0.1mm')
plt.legend(loc="lower left")

plt.show()

Un paio di considerazioni

1) LSTM e' molto sensibile all'overfitting del modello.Nel caso specifico valori di 5-6 epochs sono gia' ottimali per la convergenza del modello

2) se il modello non converge la predizione magari ha una forma e' corretta ma e' shiftata (tipo il grafico sottostante

3) si deve giocare sul valore della finestra sia come units che come time_steps (oltre che sulle epochs) per ottenere un modello ottimale

import numpy as np

import pandas as pd

import tensorflow as tf

from tensorflow.keras.models import Sequential

from tensorflow.keras.layers import LSTM, Dense, Bidirectional

from sklearn.preprocessing import MinMaxScaler

import matplotlib.pyplot as plt

data = pd.read_csv('est2.csv')

data = data['Est'].values.reshape(-1, 1)

scaler = MinMaxScaler(feature_range=(0, 1))

data_scaled = scaler.fit_transform(data)

def create_dataset(data, time_step=1):

X, y = [], []

for i in range(len(data) - time_step):

X.append(data[i:i + time_step, 0])

y.append(data[i + time_step, 0])

return np.array(X), np.array(y)

time_step = 200

X, y = create_dataset(data_scaled, time_step)

X = X.reshape(X.shape[0], X.shape[1], 1)

train_size = int(len(X) * 0.8)

X_train, X_test = X[:train_size], X[train_size:]

y_train, y_test = y[:train_size], y[train_size:]

model = Sequential()

model.add(Bidirectional(LSTM(units=200, return_sequences=True), input_shape=(time_step, 1)))

model.add(Bidirectional(LSTM(units=200)))

model.add(Dense(units=1))

model.compile(optimizer='adam', loss='mean_squared_error')

model.fit(X_train, y_train, epochs=5, batch_size=32)

predictions = model.predict(X_test)

predictions_rescaled = scaler.inverse_transform(predictions)

y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, 1))

plt.figure(figsize=(10,6))

plt.plot(y_test_rescaled, label='Actual Data')

plt.plot(predictions_rescaled, label='Predicted Data')

plt.title('Bidirectional LSTM Time Series Forecasting')

plt.xlabel('Time')

plt.ylabel('Value')

plt.legend()

plt.show()

forecast_steps = 1000 # Number of time steps to predict

input_sequence = X_test[-1]

forecast = []

for _ in range(forecast_steps):

predicted_value = model.predict(input_sequence.reshape(1, time_step, 1))

forecast.append(predicted_value[0, 0])

input_sequence = np.append(input_sequence[1:], predicted_value, axis=0)

forecast_rescaled = scaler.inverse_transform(np.array(forecast).reshape(-1, 1))

predictions_rescaled = scaler.inverse_transform(predictions)

y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, 1))

# Plot the results

plt.figure(figsize=(10,6))

plt.plot(forecast_rescaled, label='Forecasted Data')

plt.title('Bidirectional LSTM Time Series Forecasting')

plt.xlabel('Time')

plt.ylabel('Value')

plt.legend()

plt.show()

risultati sul set di validazione

questi sono i dati di forecast

Ho provato a chiedere a ChatGPT che mi ha risposto con questo codice

import numpy as np

import pandas as pd

import tensorflow as tf

from tensorflow.keras.models import Model

from tensorflow.keras.layers import LSTM, Dense, Input, Bidirectional

from sklearn.preprocessing import MinMaxScaler

from sklearn.model_selection import train_test_split

import matplotlib.pyplot as plt

def load_data(file_path):

df = pd.read_csv(file_path)

data = df['Est'].values

scaler = MinMaxScaler(feature_range=(0, 1))

data_scaled = scaler.fit_transform(data.reshape(-1, 1))

return data_scaled, scaler

def create_dataset(data, time_step=200):

X, y = [], []

for i in range(len(data) - time_step):

X.append(data[i:i + time_step, 0])

y.append(data[i + time_step, 0])

X = np.array(X)

y = np.array(y)

X = np.reshape(X, (X.shape[0], X.shape[1], 1))

return X, y

def create_bilateral_lstm_model(input_shape):

inputs = Input(shape=input_shape)

forward_lstm = LSTM(200, return_sequences=False)(inputs)

backward_lstm = LSTM(200, return_sequences=False, go_backwards=True)(inputs)

merged = tf.keras.layers.concatenate([forward_lstm, backward_lstm])

output = Dense(1)(merged)

model = Model(inputs=inputs, outputs=output)

return model

def train_and_forecast(file_path, time_step=200, forecast_days=200):

data_scaled, scaler = load_data(file_path)

X, y = create_dataset(data_scaled, time_step)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

model = create_bilateral_lstm_model(input_shape=(X_train.shape[1], 1))

model.compile(optimizer='adam', loss='mean_squared_error')

model.fit(X_train, y_train, epochs=5, batch_size=32, validation_data=(X_test, y_test))

y_pred = model.predict(X_test)

y_pred_rescaled = scaler.inverse_transform(y_pred)

y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, 1))

plt.figure(figsize=(12, 6))

plt.plot(y_test_rescaled, label="True Data")

plt.plot(y_pred_rescaled, label="Predicted Data")

plt.title("True vs Predicted Data (Test Set)")

plt.legend()

plt.show()

last_sequence = data_scaled[-time_step:].reshape(1, time_step, 1)

future_predictions = []

for _ in range(forecast_days):

forecast = model.predict(last_sequence)

future_predictions.append(forecast[0, 0])

last_sequence = np.append(last_sequence[:, 1:, :], forecast.reshape(1, 1, 1), axis=1)

future_predictions_rescaled = scaler.inverse_transform(np.array(future_predictions).reshape(-1, 1))

plt.figure(figsize=(12, 6))

plt.plot(range(len(data_scaled)), scaler.inverse_transform(data_scaled), label="Historical Data")

plt.plot(range(len(data_scaled), len(data_scaled) + forecast_days), future_predictions_rescaled, label="Forecasted Data", linestyle="--")

plt.title(f"Forecasted Data for Next {forecast_days} Steps")

plt.legend()

plt.show()

return future_predictions_rescaled

file_path = 'est2.csv'

forecast = train_and_forecast(file_path, time_step=200, forecast_days=200)

print("Forecasted future values:")

print(forecast)

ingrandimento dei soli dati di forevast

Visto che il modello converge adesso vediamo come si comporta quando si ha un evento inatteso (ovvero l'innescarsi del movimento di versante)..ho tagliato i dati sul primo movimento che corrisponda circa con la coordinata di ascissa 2000

l'addestramento del modello non e' male (anche se non perfetto)

il problema e' che quando entra in modalita' forecast sbaglia completamente indicando un movimento esattamente opposto a quello della frana

Con GRU le cose non sono andate molto bene

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

from sklearn.preprocessing import MinMaxScaler, StandardScaler

import warnings

from scipy import stats

import tensorflow as tf

from tensorflow import keras

from tensorflow.keras import Sequential, layers, callbacks

from tensorflow.keras.layers import Dense, LSTM, Dropout, GRU, Bidirectional

tf.random.set_seed(1234)

print("Num GPUs Available: ", len(tf.config.list_physical_devices('GPU')))

df = pd.read_csv("est2.csv", parse_dates = ["Date"])

# Define a function to draw time_series plot

def timeseries (x_axis, y_axis, x_label):

plt.figure(figsize = (10, 6))

plt.plot(x_axis, y_axis, color ="black")

plt.xlabel(x_label, {'fontsize': 12})

plt.ylabel('Est', {'fontsize': 12})

plt.show()

print(df)

timeseries(df.index, df['Est'], 'Time (day)')

df = df.drop('Date', axis = 1)

train_size = int(len(df)*0.8)

train_data = df.iloc[:train_size]

test_data = df.iloc[train_size:]

scaler = MinMaxScaler().fit(train_data)

train_scaled = scaler.transform(train_data)

test_scaled = scaler.transform(test_data)

def create_dataset (X, look_back = 1):

Xs, ys = [], []

for i in range(len(X)-look_back):

v = X[i:i+look_back]

Xs.append(v)

ys.append(X[i+look_back])

return np.array(Xs), np.array(ys)

LOOK_BACK = 300

X_train, y_train = create_dataset(train_scaled,LOOK_BACK)

X_test, y_test = create_dataset(test_scaled,LOOK_BACK)

# Print data shape

print("X_train.shape: ", X_train.shape)

print("y_train.shape: ", y_train.shape)

print("X_test.shape: ", X_test.shape)

print("y_test.shape: ", y_test.shape)

# Create GRU model

def create_gru(units):

model = Sequential()

# Input layer

model.add(GRU (units = units, return_sequences = True,

input_shape = [X_train.shape[1], X_train.shape[2]]))

model.add(Dropout(0.2))

# Hidden layer

model.add(GRU(units = units))

model.add(Dropout(0.2))

model.add(Dense(units = 1))

#Compile model

model.compile(optimizer='adam',loss='mse')

return model

model_gru = create_gru(300)

def fit_model(model):

early_stop = keras.callbacks.EarlyStopping(monitor = 'val_loss',patience = 10)

history = model.fit(X_train, y_train, epochs = 5,

validation_split = 0.2,

batch_size = 32, shuffle = False,

callbacks = [early_stop])

return history

history_gru = fit_model(model_gru)

y_test = scaler.inverse_transform(y_test)

y_train = scaler.inverse_transform(y_train)

def plot_loss (history, model_name):

plt.figure(figsize = (10, 6))

plt.plot(history.history['loss'])

plt.plot(history.history['val_loss'])

plt.title('Model Train vs Validation Loss for ' + model_name)

plt.ylabel('Loss')

plt.xlabel('epoch')

plt.legend(['Train loss', 'Validation loss'], loc='upper right')

plt.show()

plot_loss (history_gru, 'GRU')

# Make prediction

def prediction(model):

prediction = model.predict(X_test)

prediction = scaler.inverse_transform(prediction)

return prediction

prediction_gru = prediction(model_gru)

# Plot test data vs prediction

def plot_future(prediction, model_name, y_test):

plt.figure(figsize=(10, 6))

range_future = len(prediction)

plt.plot(np.arange(range_future), np.array(y_test), label='Test data')

plt.plot(np.arange(range_future), np.array(prediction),label='Prediction')

plt.title('Test data vs prediction for ' + model_name)

plt.legend(loc='upper left')

plt.xlabel('Time (day)')

plt.ylabel('Est')

plt.show()

plot_future(prediction_gru, 'GRU', y_test)

# bilstm 5 epoch

#model_bilstm.save("bilstm.keras")

model_gru.save("gru.keras")