0
$\begingroup$

I am training LSTM to predict tree sap flow (30 min interval). I am keeping getting negative predictions, although the most of the predictions look good. I tried to use different scalers for preprocessing but this problem consists. The full code is listed below. Simplified data is also attached. Could someone who has ever encountered this problem please help or share some experiences? Thanks! enter image description here

Data is here: https://drive.google.com/file/d/1tPACc4TAU7iYkghd4UOZFKVnfQQDxpq_/view?usp=sharing

import torch
import torch.nn as nn
import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import RobustScaler
from sklearn.model_selection import train_test_split
import csv
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt

# Define the file path for saving the CSV
csv_file_path = "best_rmse_per_clone.csv"

# Define the LSTM model
class LSTM(nn.Module):
    def __init__(self, input_size=8, hidden_layer_size=200, output_size=1):  # Change input size to 8
        super().__init__()
        self.hidden_layer_size = hidden_layer_size

        self.lstm = nn.LSTM(input_size, hidden_layer_size)

        self.linear = nn.Linear(hidden_layer_size, output_size)

        self.hidden_cell = (torch.zeros(1,1,self.hidden_layer_size),
                            torch.zeros(1,1,self.hidden_layer_size))

    def forward(self, input_seq):
        lstm_out, self.hidden_cell = self.lstm(input_seq.view(len(input_seq) ,1, -1), self.hidden_cell)
        predictions = self.linear(lstm_out.view(len(input_seq), -1))
        return predictions

# Load the data
data_origin = pd.read_csv("Sapflow_Cleaned_day153_280_simplified.csv")

# Drop NaN values and select relevant columns
selected_columns = ['Date', 'VPD', 'PAR_Den_Avg']  # Remove the clone column from selected columns
subset_data = data_origin[selected_columns].copy()  # Make a copy of the DataFrame

# Define a dictionary to store the best RMSE for each clone
best_rmse_dict = {}

# Iterate over all clones from column index 5 to the end of the dataset
for clone_index in range(5, len(data_origin.columns)):
    clone = data_origin.columns[clone_index]  # Get the clone name from column index
    clone_data = data_origin[['Date', 'VPD', 'PAR_Den_Avg', clone]].copy()  # Make a copy of the DataFrame
    
    # Drop NaN values for the current clone
    clone_data.dropna(inplace=True)
    clone_data = clone_data[clone_data[clone].values > 0.005]

    # Rename the columns for clarity
    clone_data.columns = ['Date', 'VPD', 'PAR_Den_Avg', clone]

    # Prepare data for training
    data = clone_data 
    data['Date'] = pd.to_datetime(data['Date'])

    # Extract hour from the date
    data['Hour'] = data['Date'].dt.hour
    data['Date'] = data['Date'].dt.date

    # Convert date to ordinal
    data['Date'] = data['Date'].apply(lambda x: x.toordinal())

    # Convert hour to one-hot encoding
    data = pd.get_dummies(data, columns=['Hour'], prefix='Hour')

    # Normalize the other columns
    scaler = StandardScaler()
    data['Date'] = scaler.fit_transform(data['Date'].values.reshape(-1,1))
    data['PAR_Den_Avg'] = scaler.fit_transform(data['PAR_Den_Avg'].values.reshape(-1,1))
    data['VPD'] = scaler.fit_transform(data['VPD'].values.reshape(-1,1))

    # Define the input and target variables
    X = data[['Date', 'VPD', 'PAR_Den_Avg'] + [col for col in data.columns if col.startswith('Hour_')]].values
    y = data[clone].values

    # Split the data into training and testing sets
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1, shuffle=False)

    # Convert the data into PyTorch tensors and move them to GPU
    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

    X_train_tensor = torch.tensor(X_train).float().to(device)
    y_train_tensor = torch.tensor(y_train).float().to(device)

    X_test_tensor = torch.tensor(X_test).float().to(device)
    y_test_tensor = torch.tensor(y_test).float().to(device)

    # Define the model, loss function, and optimizer and move the model to GPU
    model = LSTM(input_size=X.shape[1], hidden_layer_size=200, output_size=1).to(device)
    loss_function = nn.L1Loss()
    optimizer = torch.optim.Adam(model.parameters(), lr=0.001)

    # Train the model multiple times and select the best fit model based on RMSE
    best_rmse = float('inf')
    best_model = None

    num_trainings = 5  # Number of times to train the model
    for _ in range(num_trainings):

        # Train the model
        epochs = 40
        for i in range(epochs):
            for seq, labels in zip(X_train_tensor, y_train_tensor):
                optimizer.zero_grad()
                model.hidden_cell = (torch.zeros(1, 1, model.hidden_layer_size).to(device),
                                torch.zeros(1, 1, model.hidden_layer_size).to(device))

                y_pred = model(seq.unsqueeze(0))  # Add additional dimension

                single_loss = loss_function(y_pred, labels.unsqueeze(0))  # Adjust label shape
                single_loss.backward()
                optimizer.step()

            if i % 2 == 1:
                print(f'epoch: {i:3} loss: {single_loss.item():10.8f}')

        # Evaluate the model
        with torch.no_grad():
            test_seq = X_test_tensor.unsqueeze(1)  # Add an additional dimension
            preds = model(test_seq)
            loss = loss_function(preds, y_test_tensor.unsqueeze(1))  # Adjust label shape
            print(f'Test loss: {loss.item()}')

        # Convert the predictions back to original scale
        preds = preds.cpu().numpy()

        # Calculate RMSE
        rmse = mean_squared_error(y_test, preds, squared=False)
        print(f'RMSE: {rmse}')

        # Check if current model is the best fit
        if rmse < best_rmse:
            best_rmse = rmse
            best_model = model
    
    # Store the best RMSE for this clone in the dictionary
    best_rmse_dict[clone] = best_rmse

    # Use the best model for predictions
    with torch.no_grad():
        test_seq = X_test_tensor.unsqueeze(1)  # Add an additional dimension
        preds = best_model(test_seq)
        loss = loss_function(preds, y_test_tensor.unsqueeze(1))  # Adjust label shape
        print(f'Best model test loss: {loss.item()}')

    # Convert the predictions back to original scale
    preds = preds.cpu().numpy()

    # Plot the results
    plt.figure(figsize=(12, 6))
    plt.plot(preds, label='Predicted')
    plt.plot(y_test, label='True')
    plt.legend()
    plt.show()

    # Calculate RMSE using the best model
    rmse = mean_squared_error(y_test, preds, squared=False)
    print(f'Best model RMSE: {rmse}')

    # Save the best model state
    torch.save(best_model.state_dict(), f'D:/CNN SAP FLOW/model_{clone}_model_state_dict.pt')
    
    # Write clone and best RMSE to the CSV file
    with open(csv_file_path, mode='w', newline='') as file:
        writer = csv.writer(file)
        writer.writerow(['Clone', 'Best_RMSE'])  # Write header
        for clone, rmse in best_rmse_dict.items():
            writer.writerow([clone, rmse])
$\endgroup$

1 Answer 1

0
$\begingroup$

I solved by using BiLSTM and reducing hidden_layer_size from 100 to 60. Now it worked well as expected. enter image description here

$\endgroup$
4
  • 1
    $\begingroup$ Are you sure the bidirectional lstm isnt cheating by looking ahead? $\endgroup$
    – Karl
    Mar 30 at 17:16
  • $\begingroup$ That's a good question. The testing dataset (which is labeled as "True" here) has never been encountered during the training process. In my scenario, I aim to fill in the missing time series data to ensure a continuous "measurement" of water use. Given that trees have a growing period, my time scale is restricted from June to October. The predicted data covers only a short duration within the growing season, typically spanning just a few days. $\endgroup$
    – Jiaj
    Mar 31 at 15:50
  • $\begingroup$ BTW, consider tree's water use is driven by VPD (vapor pressure deficit) and PAR (Photosynthetic active radiation), I included them as predictors. So I think that will help to limit the "cheating", maybe? $\endgroup$
    – Jiaj
    Mar 31 at 16:02
  • $\begingroup$ It's not a matter if train/test split. If you are passing an input of multiple timesteps to a bidirectional lstm, it can use information from future timesteps to inform predictions of past timesteps. $\endgroup$
    – Karl
    Mar 31 at 17:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.