LSTM.py

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import time
import torch
from torch import nn
from torch.utils.data import DataLoader, TensorDataset
import os
import warnings
warnings.filterwarnings(
    "ignore",
    message="Converting mask without torch.bool dtype to bool; this will negatively affect performance"
)

# Setting the working directory to the current scripts location
script_dir = os.path.dirname(os.path.abspath(__file__))
os.chdir(script_dir)

# Data loading & preprocessing
def retrieve_data(loc = "./Data/"):
    """
    Iterates through four csv files containing Wind power Forcasting data and combines
    them into a single Pandas dataframe. Default location for data is the current
    directory, "Data" folder.
    """
    file_list = ["Location1.csv",
                 "Location2.csv",
                 "Location3.csv",
                 "Location4.csv"]

    data_all = []
    expected_rows = 0

    for file_id in file_list:
        file_loc = loc + file_id
        temp = pd.read_csv(file_loc)
        expected_rows += temp.shape[0]
        temp = temp.to_numpy()
        data_all.append(temp)

    data = np.concatenate(data_all, axis=0)
    print(f"expected rows: {expected_rows}")
    print(f"retrieved rows: {data.shape[0]}")

    df = pd.DataFrame(
        data,
        columns=[
            "Time","temperature_2m","relativehumidity_2m","dewpoint_2m",
            "windspeed_10m","windspeed_100m","winddirection_10m",
            "winddirection_100m","windgusts_10m","Power"
        ],
    )
    return df

def create_df(df):
    """
    Taking the DataFrame produced from retrieve_data() and modifying it for use in
    wind power forecasting. The four wind turbine locations get one-hot encoded, and
    temporal data gets converted to cyclic representation. We also added a month identifier.
    """
    df["Location"] = (
        ["LocA"] * 43800 + ["LocB"] * 43800 + ["LocC"] * 43800 + ["LocD"] * 43800
    )

    df["Time"] = pd.to_datetime(df["Time"])
    float_cols = [
        "temperature_2m","relativehumidity_2m","dewpoint_2m",
        "windspeed_10m","windspeed_100m","winddirection_10m",
        "winddirection_100m","windgusts_10m","Power"
    ]
    for c in float_cols:
        df[c] = pd.to_numeric(df[c], errors="coerce")


    df.loc[df["Location"]=="LocA","Location"]=0
    df.loc[df["Location"]=="LocB","Location"]=1
    df.loc[df["Location"]=="LocC","Location"]=2
    df.loc[df["Location"]=="LocD","Location"]=3
    df["Location"] = df["Location"].astype(int)
    oh = pd.get_dummies(df["Location"])
    df["LocA"], df["LocB"], df["LocC"], df["LocD"] = oh[0], oh[1], oh[2], oh[3]


    df["Hour"]  = df["Time"].dt.hour
    df["Day"]   = df["Time"].dt.day
    df["Month"] = df["Time"].dt.month
    df["Year"]  = df["Time"].dt.year
    df["Hour"] = (df["Hour"]/24.0) * 2*np.pi
    df["hour_sin"]  = np.sin(df["Hour"]);  df["hour_cos"]  = np.cos(df["Hour"])
    df["day_sin"]   = np.sin(2*np.pi*df["Day"]/31.0);  df["day_cos"] = np.cos(2*np.pi*df["Day"]/31.0)
    df["month_sin"] = np.sin(2*np.pi*df["Month"]/12.0); df["month_cos"] = np.cos(2*np.pi*df["Month"]/12.0)
    df.drop(columns=["Hour","Time","Day","Month"], inplace=True)
    print(f"Original Shape:{df.shape}")
    return df

features = [
    'temperature_2m','relativehumidity_2m','dewpoint_2m',
    'windspeed_10m','windspeed_100m','winddirection_10m',
    'winddirection_100m','windgusts_10m',
    'LocA','LocB','LocC','LocD',
    'hour_sin','hour_cos','day_sin','day_cos','month_sin','month_cos',
    'Power'
]
scale_features = [
    'temperature_2m','relativehumidity_2m','dewpoint_2m',
    'windspeed_10m','windspeed_100m','winddirection_10m',
    'winddirection_100m','windgusts_10m',
    'LocA','LocB','LocC','LocD',
    'hour_sin','hour_cos','day_sin','day_cos','month_sin','month_cos'
]

target_name = 'Power'
def split_scale_per_location(df, train_frac=0.7, val_frac=0.15):
    trains, vals, tests = [], [], []
    for loc_col in ["LocA","LocB","LocC","LocD"]:
        g = df[df[loc_col]==1].reset_index(drop=True)
        n = len(g)
        i_train = int(n*train_frac)
        i_val   = int(n*(train_frac+val_frac))
        trains.append(g.iloc[:i_train])
        vals.append(g.iloc[i_train:i_val])
        tests.append(g.iloc[i_val:])
    df_train = pd.concat(trains, ignore_index=True)
    df_val   = pd.concat(vals, ignore_index=True)
    df_test  = pd.concat(tests, ignore_index=True)

    scaler = MinMaxScaler()
    df_train[scale_features] = scaler.fit_transform(df_train[scale_features])
    df_val[scale_features]   = scaler.transform(df_val[scale_features])
    df_test[scale_features]  = scaler.transform(df_test[scale_features])
    return df_train, df_val, df_test, scaler

def make_windows(ndarr, input_len=24, output_len=6, target_idx=None):
    X, y = [], []
    for i in range(len(ndarr)-input_len-output_len+1):
        X.append(ndarr[i:i+input_len])
        y.append(ndarr[i+input_len:i+input_len+output_len, target_idx])
    return np.array(X), np.array(y)

def build_windows_all_locs(df, input_len=24, output_len=6):
    tgt_idx = features.index(target_name)
    Xs, ys = [], []
    for loc_col in ["LocA","LocB","LocC","LocD"]:
        g = df[df[loc_col]==1].reset_index(drop=True)
        Xg, yg = make_windows(g[features].to_numpy(), input_len, output_len, tgt_idx)
        Xs.append(Xg); ys.append(yg)
    return np.concatenate(Xs), np.concatenate(ys)

def make_loaders(Xtr,Ytr,Xva,Yva,Xte,Yte,batch_size=64,num_workers=0):
    def to_ds(X,y): return TensorDataset(torch.from_numpy(X).float(), torch.from_numpy(y).float())
    tr, va, te = to_ds(Xtr,Ytr), to_ds(Xva,Yva), to_ds(Xte,Yte)
    train_loader = DataLoader(tr, batch_size=batch_size, shuffle=True, drop_last=True,  num_workers=num_workers)
    val_loader   = DataLoader(va, batch_size=batch_size, shuffle=False, drop_last=False, num_workers=num_workers)
    test_loader  = DataLoader(te, batch_size=batch_size, shuffle=False, drop_last=False, num_workers=num_workers)
    return train_loader, val_loader, test_loader


# LSTM Encoder–Decoder
class LSTMEncoderDecoder(nn.Module):
    """
    Autoregressive sequence-to-sequence LSTM:
      - Encoder reads [B, T_in, D] and returns hidden/cell
      - Decoder generates output steps
      - Optional teacher forcing during training
    """
    def __init__(self, input_size, hidden_size=128, num_layers=2, output_len=6, dropout=0.1):
        super().__init__()
        self.encoder = nn.LSTM(input_size, hidden_size, num_layers,
                               batch_first=True, dropout=dropout)
        self.decoder = nn.LSTM(1, hidden_size, num_layers,
                               batch_first=True, dropout=dropout)
        self.fc = nn.Linear(hidden_size, 1)
        self.output_len = output_len

    def forward(self, x, tgt_seq=None, teacher_forcing_ratio=0.0):

        B = x.size(0)
        _, (h, c) = self.encoder(x)
        last_power = x[:, -1:, -1:]
        dec_in = last_power
        outputs = []
        use_tf = (tgt_seq is not None) and (teacher_forcing_ratio > 0)
        if use_tf and tgt_seq.dim()==2:
            tgt_seq = tgt_seq.unsqueeze(-1)

        for t in range(self.output_len):
            out, (h, c) = self.decoder(dec_in, (h, c))
            pred = self.fc(out)
            outputs.append(pred)

            if use_tf and (np.random.rand() < teacher_forcing_ratio):
                next_in = tgt_seq[:, t:t+1, :]  # ground-truth step t
            else:
                next_in = pred
            dec_in = next_in

        return torch.cat(outputs, dim=1).squeeze(-1)  # [B, T_out]


# Training / Evaluation utils

def calculate_metrics(y_true, y_pred):
    mae = mean_absolute_error(y_true.flatten(), y_pred.flatten())
    mse = mean_squared_error(y_true.flatten(), y_pred.flatten())
    rmse = np.sqrt(mse)
    r2 = r2_score(y_true.flatten(), y_pred.flatten())
    return {'MAE': mae, 'MSE': mse, 'RMSE': rmse, 'R2': r2}

class EarlyStopping:
    def __init__(self, patience=10, min_delta=1e-6, restore_best_weights=True):
        self.patience = patience
        self.min_delta = min_delta
        self.restore_best_weights = restore_best_weights
        self.best_loss = None
        self.counter = 0
        self.best_state = None

    def step(self, val_loss, model):
        if self.best_loss is None or val_loss < self.best_loss - self.min_delta:
            self.best_loss = val_loss
            self.counter = 0
            if self.restore_best_weights:
                self.best_state = {k: v.detach().cpu().clone()
                                   for k, v in model.state_dict().items()}
        else:
            self.counter += 1
        stop = self.counter >= self.patience
        if stop and self.restore_best_weights and self.best_state is not None:
            model.load_state_dict(self.best_state)
        return stop

def train_model(model, train_loader, val_loader, device,
                num_epochs=50, lr=1e-3, teacher_forcing_ratio=0.5, print_every=1):
    criterion = nn.HuberLoss(delta=1.0)
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)

    es = EarlyStopping(patience=10, min_delta=1e-6)

    train_losses, val_losses = [], []
    for epoch in range(1, num_epochs+1):

        model.train()
        tr_loss = 0.0
        for xb, yb in train_loader:
            xb, yb = xb.to(device), yb.to(device)
            optimizer.zero_grad(set_to_none=True)
            preds = model(xb, tgt_seq=yb, teacher_forcing_ratio=teacher_forcing_ratio)
            loss = criterion(preds, yb)
            loss.backward()
            nn.utils.clip_grad_norm_(model.parameters(), 1.0)
            optimizer.step()
            tr_loss += loss.item()
        tr_loss /= len(train_loader)
        train_losses.append(tr_loss)


        model.eval()
        va_loss = 0.0
        with torch.no_grad():
            for xb, yb in val_loader:
                xb, yb = xb.to(device), yb.to(device)
                preds = model(xb, tgt_seq=None, teacher_forcing_ratio=0.0)
                va_loss += criterion(preds, yb).item()
        va_loss /= len(val_loader)
        val_losses.append(va_loss)

        if epoch % print_every == 0:
            print(f"Epoch {epoch}/{num_epochs} - Train Loss: {tr_loss:.6f} | Val Loss: {va_loss:.6f}")

        if es.step(va_loss, model):
            print(f"Early stopping at epoch {epoch}.")
            break

    return train_losses, val_losses


# Main
if __name__ == "__main__":
    torch.manual_seed(42)
    np.random.seed(42)

    # edit if required
    data_loc = "./Data/"

    pred_horizon_list = [2,4,6,8,10,12]
    metrics_list = []
    time_list = []
    
    
    for pred_horizon in pred_horizon_list:
    
        start = time.time()
        
        # Hyperparameters
        input_len = 24
        output_len = pred_horizon
        batch_size = 64
        learning_rate = 1e-3
        epochs = 200
        teacher_forcing = 0.5

        device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
        print("Using device:", device)

        print("Loading & preparing data")
        raw = retrieve_data(loc = data_loc)
        df = create_df(raw)
        df = df.drop(columns=["Location"])

        df_train, df_val, df_test, scaler = split_scale_per_location(df)

        Xtr, Ytr = build_windows_all_locs(df_train, input_len, output_len)
        Xva, Yva = build_windows_all_locs(df_val,   input_len, output_len)
        Xte, Yte = build_windows_all_locs(df_test,  input_len, output_len)

        print("Shapes:")
        print("  train:", Xtr.shape, Ytr.shape)
        print("  val  :", Xva.shape, Yva.shape)
        print("  test :", Xte.shape, Yte.shape)

        train_loader, val_loader, test_loader = make_loaders(
            Xtr, Ytr, Xva, Yva, Xte, Yte, batch_size=batch_size, num_workers=0
        )

        # Model
        model = LSTMEncoderDecoder(
            input_size=Xtr.shape[-1],
            hidden_size=128,
            num_layers=2,
            output_len=output_len,
            dropout=0.1
        ).to(device)

        print("\nTraining")
        train_losses, val_losses = train_model(
            model, train_loader, val_loader, device,
            num_epochs=epochs, lr=learning_rate, teacher_forcing_ratio=teacher_forcing
        )

        print("\nEvaluating on test set")
        model.eval()
        test_preds, test_tgts = [], []
        criterion = nn.HuberLoss(delta=1.0)
        test_loss = 0.0
        with torch.no_grad():
            for xb, yb in test_loader:
                xb, yb = xb.to(device), yb.to(device)
                preds = model(xb, tgt_seq=None, teacher_forcing_ratio=0.0)
                test_preds.append(preds.cpu().numpy())
                test_tgts.append(yb.cpu().numpy())
                test_loss += criterion(preds, yb).item()
        test_loss /= len(test_loader)
        test_preds = np.concatenate(test_preds, axis=0)
        test_tgts  = np.concatenate(test_tgts,  axis=0)

        metrics = calculate_metrics(test_tgts, test_preds)
        print("\nTest Set Performance:")
        for k,v in metrics.items():
            print(f"{k}: {v:.6f}")

        end = time.time()
        time_taken = end-start
        
        metrics_list.append(metrics)
        time_list.append(time_taken)
        
        
        print(f"Execution time: {end - start:.4f} seconds")

    # Final results for all time horizons
    df = pd.DataFrame({'pred_length': pred_horizon_list,
                          'metrics' : metrics_list,
                            'time' : time_list})
                            
    df.to_csv('LSTM_results.csv')

    # Train/Val Loss Plot
    plt.figure(figsize=(12,4))
    plt.plot(train_losses, label="Train Loss")
    plt.plot(val_losses,   label="Val Loss")
    plt.xlabel("Epoch"); plt.ylabel("Huber Loss")
    plt.title("Training & Validation Loss")
    plt.grid(True); plt.legend(); plt.tight_layout()
    plt.savefig("lstm_train_val_loss.png")

    # Actual vs Predicted Plot
    plt.figure(figsize=(5.5,5.5))
    plt.scatter(test_tgts.flatten(), test_preds.flatten(), alpha=0.35)
    lo, hi = test_tgts.min(), test_tgts.max()
    plt.plot([lo,hi],[lo,hi],'r--')
    plt.xlabel("Actual"); plt.ylabel("Predicted")
    plt.title(f"Actual vs Predicted (R2={metrics['R2']:.3f})")
    plt.grid(True); plt.tight_layout()
    plt.savefig("lstm_actual_vs_pred.png")