real

README.md

@@ -1,28 +1,5 @@
 # midterm ml&dm1
 
-## not cool
+## real
 
-currently this shit is a blatant copy. changing soon.
-
-## how to
-
-load this thing on google colab dont run on device its so cooked
-
-how to git? do not `push --force`. prefer pull requests. its like github
-
-## roadmap
-
-we will plan to do this in multiple different ways using scikit-learn and pytorch. tensorflow is ancient.
-
-> [!NOTE]
-> ~phase out tensorflow. dont use it.~
->
-> its 2 different platforms for 2 different types of people. use whatever ig
-
-maybe also fit more weather related stats? this dataset is taken from nasa so sources are there
-
-## collab
-
-seng tells me this project can be planned without code. basically someone find da waes and other people implements them
-
-thing is you click [here at scikit learn](https://scikit-learn.org/stable/supervised_learning.html) and tell me how many ways there are you can use to solve a problem. holy fucking shit. too many.
+not real

data_loader.py

@@ -32,9 +32,10 @@ def load_weather_data(file_path: str) -> pd.DataFrame:
 
     # Convert YEAR, MO, DY to datetime index
     if all(col in data.columns for col in ("YEAR", "MO", "DY")):
-        data["date"] = pd.to_datetime(
-            {k: data[v] for v, k in zip(["YEAR", "MO", "DY"], ["year", "month", "day"])}
-        )
+        dt_conv_table = {
+            k: data[v] for v, k in zip(["YEAR", "MO", "DY"], ["year", "month", "day"])
+        }
+        data["date"] = pd.to_datetime(dt_conv_table)
         data.drop(columns=["YEAR", "MO", "DY"], inplace=True)
         data.set_index("date", inplace=True)
 

dataset_update.py

@@ -6,7 +6,7 @@ data_file = "hanoi_weather_data.csv"
 
 # Define parameters
 latitude, longitude = 21.04815477090267, 105.80123159165096  # lmao
-start_date = "1985-01-01"
+start_date = "1995-01-01"
 end_date = datetime.now().strftime("%Y-%m-%d")  # Today's date
 
 # Define desired parameters
@@ -19,7 +19,10 @@ end_date = datetime.now().strftime("%Y-%m-%d")  # Today's date
 # WD10M           MERRA-2 Wind Direction at 10 Meters (Degrees)
 # more at https://power.larc.nasa.gov/data-access-viewer/
 weather_params = [
+    "T2M",
     "QV2M",
+    "RH2M",
+    "PS",
 ]
 
 params = {

dead/dead/pytorch_port.py

@@ -1,163 +0,0 @@
-from sklearn.model_selection import train_test_split
-import torch
-
-device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-
-from sklearn.preprocessing import MinMaxScaler
-
-import numpy as np
-from load_data import load_this_data
-
-a = load_this_data()
-
-# time series forecasting
-
-
-scaler = MinMaxScaler()
-print(a.columns)
-scaled_data = scaler.fit_transform(a)
-
-
-def create_sequences(data, seq_length):
-    "see this? return[0] is the big x and return[1] is the big y. and h(x) = y. and h(x) is the model."
-    xs = []
-    ys = []
-    for i in range(len(data) - seq_length):
-        x = data[i : (i + seq_length)]
-        y = data[i + seq_length]
-        xs.append(x)
-        ys.append(y)
-    return np.array(xs), np.array(ys)
-
-
-seq_length = 20
-
-X, y = create_sequences(scaled_data, seq_length)
-
-X_train, X_test, y_train, y_test = train_test_split(
-    X, y, test_size=0.2, random_state=42
-)
-from torch.utils.data import TensorDataset, DataLoader
-
-X_train_tensor = torch.tensor(X_train, dtype=torch.float32)
-y_train_tensor = torch.tensor(y_train, dtype=torch.float32)
-
-train_dataset = TensorDataset(X_train_tensor, y_train_tensor)
-train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)
-
-
-import torch.nn as nn
-
-
-class LSTMModel(nn.Module):
-    def __init__(self, input_size, hidden_size, num_layers, output_size):
-        super(LSTMModel, self).__init__()
-        self.hidden_size = hidden_size
-        self.num_layers = num_layers
-        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
-        self.fc = nn.Linear(hidden_size, output_size)
-
-    def forward(self, x):
-        h0 = torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(device)
-        c0 = torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(device)
-
-        out, _ = self.lstm(x, (h0, c0))
-        out = self.fc(out[:, -1, :])  # Take the output from the last time step
-        return out
-
-
-model = LSTMModel(
-    input_size=X_train.shape[2],
-    hidden_size=64,
-    num_layers=2,
-    output_size=y_train.shape[1],
-)
-model.to(device)
-criterion = nn.MSELoss()
-optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
-
-train_losses = []
-val_losses = []
-
-
-for epoch in range(100):
-    model.train()
-    running_train_loss = 0.0
-
-    # inputs, targets
-    for batch_X, batch_y in train_loader:
-        batch_X = batch_X.to(device)
-        batch_y = batch_y.to(device)
-
-        outputs = model(batch_X)
-        loss = criterion(outputs, batch_y)
-
-        optimizer.zero_grad()
-        loss.backward()
-        optimizer.step()
-
-        running_train_loss += loss.item() * batch_X.size(0)
-    if (epoch + 1) % 10 == 0:
-        print(f"Epoch [{epoch+1}/100], Loss: {loss.item():.4f}")
-
-    epoch_train_loss = running_train_loss / len(train_loader.dataset)
-    train_losses.append(epoch_train_loss)
-
-model.eval()
-
-X_test_tensor = torch.tensor(X_test, dtype=torch.float32).to(device)
-y_test_tensor = torch.tensor(y_test, dtype=torch.float32).to(device)
-val_dataset = TensorDataset(X_test_tensor, y_test_tensor)
-val_loader = DataLoader(val_dataset, batch_size=64, shuffle=False)
-
-# with torch.no_grad():
-#     running_val_loss = 0.0
-#     predictions = model(X_test_tensor)
-#     test_loss = criterion(predictions, y_test_tensor)
-#     print(f"Test Loss: {test_loss.item():.4f}")
-
-running_val_loss = 0.0
-predictions = []
-with torch.no_grad():
-    for inputs, targets in val_loader:
-        inputs, targets = inputs.to(device), targets.to(device)
-        outputs = model(inputs)
-        loss = criterion(outputs, targets)
-        running_val_loss += loss.item() * inputs.size(0)
-        predictions.append(outputs.cpu().numpy())  # Collect all batches
-
-epoch_val_loss = running_val_loss / len(val_loader.dataset)  # len(val_dataset)
-val_losses.append(epoch_val_loss)
-
-predictions = np.concatenate(predictions, axis=0)
-import matplotlib.pyplot as plt
-
-y_test_cpu = y_test_tensor.cpu().numpy()
-
-predictions_rescaled = scaler.inverse_transform(predictions)
-y_test_rescaled = scaler.inverse_transform(y_test_cpu)
-
-
-fig, (ax1, ax2) = plt.subplots(
-    1, 2, figsize=(14, 5), gridspec_kw={"width_ratios": [1, 3]}
-)
-
-# Loss curve on the left
-ax1.plot(train_losses, label="Train Loss")
-ax1.plot(val_losses, label="Val Loss")
-ax1.set_title("Loss Curve")
-ax1.set_xlabel("Epoch")
-ax1.set_ylabel("Loss")
-ax1.legend()
-ax1.grid(True)
-
-# Prediction vs Actual on the right
-ax2.plot(y_test_rescaled, label="Actual")
-ax2.plot(predictions_rescaled, label="Predicted")
-ax2.set_title("Actual vs. Predicted")
-ax2.set_xlabel("Sample")
-ax2.set_ylabel("Value")
-ax2.legend()
-
-plt.tight_layout()
-plt.show()

dead/global_config.py

@@ -1 +0,0 @@
-DATA_FILE_NAME = "hanoi_weather_data.csv"

dead/load_data.py

@@ -1,28 +0,0 @@
-"extracting the load data step from the notebook"
-
-import io
-
-# from itertools import islice
-from typing import Union
-
-def cut_headers(stream: Union[io.StringIO, io.FileIO, io.TextIOWrapper]):
-    headers = []
-    exist_startheader = False
-
-    while line := stream.readline():
-        if "-BEGIN HEADER-" in line and not exist_startheader:
-            exist_startheader = True
-            headers.append(str(line.strip()))
-        elif "-END HEADER-" in line and exist_startheader:
-            headers.append(str(line.strip()))
-            break
-        elif exist_startheader:
-            # print(f"{line = }")
-            headers.append(str(line.strip()))
-        elif not exist_startheader:
-            continue
-    if not exist_startheader:
-        # If no headers are found, reset the stream and return
-        stream.seek(0)
-        return "", stream
-    return "\n".join(headers), stream

dead/pytorch_port.py

@@ -1,163 +0,0 @@
-from sklearn.model_selection import train_test_split
-import torch
-
-device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
-
-from sklearn.preprocessing import MinMaxScaler
-
-import numpy as np
-from load_data import load_this_data
-
-a = load_this_data()
-
-# time series forecasting
-
-
-scaler = MinMaxScaler()
-print(a.columns)
-scaled_data = scaler.fit_transform(a)
-
-
-def create_sequences(data, seq_length):
-    "see this? return[0] is the big x and return[1] is the big y. and h(x) = y. and h(x) is the model."
-    xs = []
-    ys = []
-    for i in range(len(data) - seq_length):
-        x = data[i : (i + seq_length)]
-        y = data[i + seq_length]
-        xs.append(x)
-        ys.append(y)
-    return np.array(xs), np.array(ys)
-
-
-seq_length = 20
-
-X, y = create_sequences(scaled_data, seq_length)
-
-X_train, X_test, y_train, y_test = train_test_split(
-    X, y, test_size=0.2, random_state=42
-)
-from torch.utils.data import TensorDataset, DataLoader
-
-X_train_tensor = torch.tensor(X_train, dtype=torch.float32)
-y_train_tensor = torch.tensor(y_train, dtype=torch.float32)
-
-train_dataset = TensorDataset(X_train_tensor, y_train_tensor)
-train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)
-
-
-import torch.nn as nn
-
-
-class LSTMModel(nn.Module):
-    def __init__(self, input_size, hidden_size, num_layers, output_size):
-        super(LSTMModel, self).__init__()
-        self.hidden_size = hidden_size
-        self.num_layers = num_layers
-        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
-        self.fc = nn.Linear(hidden_size, output_size)
-
-    def forward(self, x):
-        h0 = torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(device)
-        c0 = torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(device)
-
-        out, _ = self.lstm(x, (h0, c0))
-        out = self.fc(out[:, -1, :])  # Take the output from the last time step
-        return out
-
-
-model = LSTMModel(
-    input_size=X_train.shape[2],
-    hidden_size=64,
-    num_layers=2,
-    output_size=y_train.shape[1],
-)
-model.to(device)
-criterion = nn.MSELoss()
-optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
-
-train_losses = []
-val_losses = []
-
-
-for epoch in range(100):
-    model.train()
-    running_train_loss = 0.0
-
-    # inputs, targets
-    for batch_X, batch_y in train_loader:
-        batch_X = batch_X.to(device)
-        batch_y = batch_y.to(device)
-
-        outputs = model(batch_X)
-        loss = criterion(outputs, batch_y)
-
-        optimizer.zero_grad()
-        loss.backward()
-        optimizer.step()
-
-        running_train_loss += loss.item() * batch_X.size(0)
-    if (epoch + 1) % 10 == 0:
-        print(f"Epoch [{epoch+1}/100], Loss: {loss.item():.4f}")
-
-    epoch_train_loss = running_train_loss / len(train_loader.dataset)
-    train_losses.append(epoch_train_loss)
-
-model.eval()
-
-X_test_tensor = torch.tensor(X_test, dtype=torch.float32).to(device)
-y_test_tensor = torch.tensor(y_test, dtype=torch.float32).to(device)
-val_dataset = TensorDataset(X_test_tensor, y_test_tensor)
-val_loader = DataLoader(val_dataset, batch_size=64, shuffle=False)
-
-# with torch.no_grad():
-#     running_val_loss = 0.0
-#     predictions = model(X_test_tensor)
-#     test_loss = criterion(predictions, y_test_tensor)
-#     print(f"Test Loss: {test_loss.item():.4f}")
-
-running_val_loss = 0.0
-predictions = []
-with torch.no_grad():
-    for inputs, targets in val_loader:
-        inputs, targets = inputs.to(device), targets.to(device)
-        outputs = model(inputs)
-        loss = criterion(outputs, targets)
-        running_val_loss += loss.item() * inputs.size(0)
-        predictions.append(outputs.cpu().numpy())  # Collect all batches
-
-epoch_val_loss = running_val_loss / len(val_loader.dataset)  # len(val_dataset)
-val_losses.append(epoch_val_loss)
-
-predictions = np.concatenate(predictions, axis=0)
-import matplotlib.pyplot as plt
-
-y_test_cpu = y_test_tensor.cpu().numpy()
-
-predictions_rescaled = scaler.inverse_transform(predictions)
-y_test_rescaled = scaler.inverse_transform(y_test_cpu)
-
-
-fig, (ax1, ax2) = plt.subplots(
-    1, 2, figsize=(14, 5), gridspec_kw={"width_ratios": [1, 3]}
-)
-
-# Loss curve on the left
-ax1.plot(train_losses, label="Train Loss")
-ax1.plot(val_losses, label="Val Loss")
-ax1.set_title("Loss Curve")
-ax1.set_xlabel("Epoch")
-ax1.set_ylabel("Loss")
-ax1.legend()
-ax1.grid(True)
-
-# Prediction vs Actual on the right
-ax2.plot(y_test_rescaled, label="Actual")
-ax2.plot(predictions_rescaled, label="Predicted")
-ax2.set_title("Actual vs. Predicted")
-ax2.set_xlabel("Sample")
-ax2.set_ylabel("Value")
-ax2.legend()
-
-plt.tight_layout()
-plt.show()

dead/sk2_port.py

@@ -1,185 +0,0 @@
-import sys
-import numpy as np
-import pandas as pd
-import matplotlib.pyplot as plt
-import seaborn as sns
-import io
-from typing import Union
-from sklearn.metrics import mean_squared_error, r2_score
-from sklearn.preprocessing import MinMaxScaler
-from sklearn.model_selection import train_test_split
-from sklearn.linear_model import LinearRegression
-from statsmodels.tsa.ar_model import AutoReg
-
-# ==== LOAD & CLEAN DATA ====
-DATA_FILE_NAME = "hanoi_weather_data.csv"
-
-def cut_headers(stream: Union[io.StringIO, io.FileIO, io.TextIOWrapper]):
-    headers = []
-    exist_startheader = False
-    while line := stream.readline():
-        if "-BEGIN HEADER-" in line and not exist_startheader:
-            exist_startheader = True
-            headers.append(str(line.strip()))
-        elif "-END HEADER-" in line and exist_startheader:
-            headers.append(str(line.strip()))
-            break
-        elif exist_startheader:
-            headers.append(str(line.strip()))
-        elif not exist_startheader:
-            continue
-    if not exist_startheader:
-        stream.seek(0)
-        return "", stream
-    return "\n".join(headers), stream
-
-with open(DATA_FILE_NAME, "r", encoding="utf-8") as file:
-    data_header, csv_file = cut_headers(file)
-    data = pd.read_csv(csv_file, na_values=-999)
-
-if all(z in data.columns for z in ("YEAR", "MO", "DY")):
-    data["date"] = pd.to_datetime(
-        {k: data[v] for v, k in zip(["YEAR", "MO", "DY"], ["year", "month", "day"])}
-    )
-    data.drop(columns=["YEAR", "MO", "DY"], inplace=True)
-    data.set_index("date", inplace=True)
-
-data.dropna(inplace=True)
-
-# ==== VISUALIZE CORRELATION ====
-plt.figure(figsize=(12, 8))
-sns.heatmap(data.corr(), annot=True, cmap="coolwarm")
-plt.title("Feature Correlation Heatmap")
-plt.show()
-
-# ==== SCALE DATA ====
-scaler = MinMaxScaler()
-scaled_data = scaler.fit_transform(data)
-
-# ==== CREATE SEQUENCES ====
-def create_sequences(data, seq_length):
-    xs, ys = [], []
-    for i in range(len(data) - seq_length):
-        x = data[i:i + seq_length].flatten()
-        y = data[i + seq_length]
-        xs.append(x)
-        ys.append(y)
-    return np.array(xs), np.array(ys)
-
-seq_length = 50
-X, y = create_sequences(scaled_data, seq_length)
-
-# ==== SPLIT ====
-date_index = data.index[seq_length:]
-X_train, X_test, y_train, y_test, date_train, date_test = train_test_split(
-    X, y, date_index, test_size=0.02, shuffle=False
-)
-
-# ==== TRAIN LINEAR REGRESSION ====
-model = LinearRegression()
-model.fit(X_train, y_train)
-predictions = model.predict(X_test)
-
-# ==== EVALUATION ====
-mse = mean_squared_error(y_test, predictions)
-r2 = r2_score(y_test, predictions)
-
-print(f"Test MSE: {mse:.4f}")
-print(f"Test R2: {r2:.4f}")
-print(f"Model Score: {model.score(X_test, y_test):.4f}")
-
-# ==== INVERSE TRANSFORM ====
-predictions_rescaled = scaler.inverse_transform(
-    predictions.reshape(-1, data.shape[1])
-)
-y_test_rescaled = scaler.inverse_transform(
-    y_test.reshape(-1, data.shape[1])
-)
-
-# ==== PLOTS ====
-
-# Actual vs Predicted with Real Dates
-plt.figure(figsize=(12, 6))
-plt.plot(date_test, y_test_rescaled[:, 0], label="Actual")
-plt.plot(date_test, predictions_rescaled[:, 0], label="Predicted")
-plt.title("Actual vs Predicted (First Feature)")
-plt.xlabel("Date")
-plt.ylabel("Original Scale")
-plt.xticks(rotation=45)
-plt.legend()
-plt.grid(True)
-plt.tight_layout()
-plt.show()
-
-# Residual Plot
-residuals = y_test[:, 0] - predictions[:, 0]
-plt.figure(figsize=(10, 4))
-plt.plot(date_test, residuals, label="Residuals")
-plt.hlines(0, 0, len(residuals), colors="red")
-plt.title("Residuals Over Time (First Feature)")
-plt.xlabel("Date")
-plt.ylabel("Error")
-plt.xticks(rotation=45)
-plt.legend()
-plt.tight_layout()
-plt.grid(True)
-plt.show()
-
-# Error Distribution
-errors = np.abs(y_test[:, 0] - predictions[:, 0])
-plt.hist(errors, bins=50)
-plt.title("Error Distribution (First Feature)")
-plt.xlabel("Absolute Error")
-plt.ylabel("Frequency")
-plt.grid(True)
-plt.show()
-
-# Scatter Plot
-plt.scatter(y_test[:, 0], predictions[:, 0], alpha=0.5)
-plt.plot([0, 1], [0, 1], "r--")
-plt.title("Predicted vs Actual (Normalized)")
-plt.xlabel("Actual")
-plt.ylabel("Predicted")
-plt.grid(True)
-plt.show()
-
-# ==== AUTOREGRESSION FORECAST ====
-target_column = data.columns[0]
-series = data[target_column]
-
-train_size = int(len(series) * 0.98)
-train_series = series[:train_size]
-test_series = series[train_size:]
-
-ar_model = AutoReg(train_series, lags=30).fit()
-ar_preds = ar_model.predict(start=len(train_series), end=len(series)-1)
-
-# Plot AR Forecast
-plt.figure(figsize=(12, 5))
-plt.plot(test_series.index, test_series, label="Actual")
-plt.plot(test_series.index, ar_preds, label="AR Forecast", linestyle="--")
-plt.title(f"AutoRegression Forecast - {target_column}")
-plt.xlabel("Date")
-plt.ylabel(target_column)
-plt.xticks(rotation=45)
-plt.legend()
-plt.grid(True)
-plt.tight_layout()
-plt.show()
-
-# ==== FUTURE FORECAST ====
-future_steps = 30
-future_index = pd.date_range(start=series.index[-1], periods=future_steps+1, freq="D")[1:]
-future_forecast = ar_model.predict(start=len(series), end=len(series) + future_steps - 1)
-
-plt.figure(figsize=(12, 5))
-plt.plot(series[-100:], label="Recent Data")
-plt.plot(future_index, future_forecast, label="Future Forecast", color="red")
-plt.title(f"{target_column} - 30 Day Forecast (AutoRegression)")
-plt.xlabel("Date")
-plt.ylabel(target_column)
-plt.xticks(rotation=45)
-plt.legend()
-plt.tight_layout()
-plt.grid(True)
-plt.show()

dead/sklearn_port.py

@@ -1,110 +0,0 @@
-import sys
-import numpy as np
-import pandas as pd
-from sklearn.metrics import r2_score
-from sklearn.preprocessing import MinMaxScaler
-from sklearn.model_selection import train_test_split
-from sklearn.linear_model import LinearRegression
-from sklearn.metrics import mean_squared_error
-import matplotlib.pyplot as plt
-# from global_config import DATA_FILE_NAME
-DATA_FILE_NAME = "hanoi_weather_data.csv"
-
-# from load_data import cut_headers
-
-import io
-
-# from itertools import islice
-from typing import Union
-
-def cut_headers(stream: Union[io.StringIO, io.FileIO, io.TextIOWrapper]):
-    headers = []
-    exist_startheader = False
-
-    while line := stream.readline():
-        if "-BEGIN HEADER-" in line and not exist_startheader:
-            exist_startheader = True
-            headers.append(str(line.strip()))
-        elif "-END HEADER-" in line and exist_startheader:
-            headers.append(str(line.strip()))
-            break
-        elif exist_startheader:
-            # print(f"{line = }")
-            headers.append(str(line.strip()))
-        elif not exist_startheader:
-            continue
-    if not exist_startheader:
-        # If no headers are found, reset the stream and return
-        stream.seek(0)
-        return "", stream
-    return "\n".join(headers), stream
-
-# 1. Load data
-
-with open(DATA_FILE_NAME, "r", encoding="utf-8") as file:
-    data_header, csv_file = cut_headers(file)
-    data = pd.read_csv(csv_file, na_values=-999)
-
-if all(z in data.columns for z in ("YEAR", "MO", "DY")):
-    data["date"] = pd.to_datetime(
-        {k: data[v] for v, k in (zip(["YEAR", "MO", "DY"], ["year", "month", "day"]))}
-    )
-    data.drop(inplace=True, columns=["YEAR", "MO", "DY"])
-    data.set_index(["date"], inplace=True)
-data.dropna(inplace=True)
-
-# 2. Normalize data
-scaler = MinMaxScaler()
-print(data.columns)
-scaled_data = scaler.fit_transform(data)
-
-
-# 3. Create sequences for supervised learning
-def create_sequences(data, seq_length):
-    xs, ys = [], []
-    for i in range(len(data) - seq_length):
-        x = data[i : i + seq_length].flatten()
-        y = data[i + seq_length]
-        xs.append(x)
-        ys.append(y)
-    return np.array(xs), np.array(ys)
-
-
-seq_length = 50
-X, y = create_sequences(scaled_data, seq_length)
-
-# 4. Split into training and testing
-X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.02)
-
-# 5. Train linear regression model
-model = LinearRegression()
-model.fit(X_train, y_train)
-
-# 6. Make predictions
-predictions = model.predict(X_test)
-
-# 7. Evaluate performance
-mse = mean_squared_error(y_test, predictions)
-r2 = r2_score(y_test, predictions)
-
-print(f"Test MSE: {mse:.4f}")
-print(f"Test r2: {r2:.4f}")
-
-print(f"Score: {model.score(X_test, y_test)}")
-# 8. Inverse transform predictions
-predictions_rescaled = scaler.inverse_transform(predictions.reshape(-1, data.shape[1]))
-y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, data.shape[1]))
-
-# 9. Plot results
-
-# Loss bar
-
-
-# Prediction vs Actual
-plt.plot(y_test_rescaled, label="Actual")
-plt.plot(predictions_rescaled, label="Predicted")
-plt.title("Actual vs. Predicted")
-plt.xlabel("Sample")
-plt.ylabel("Value")
-plt.legend()
-plt.show()

main.py

@@ -1,28 +1,58 @@
 # import pandas as pd
+import os
+from matplotlib import pyplot as plt
+import numpy as np
+from sklearn.model_selection import train_test_split  # type: ignore[import]
+
 from data_loader import load_weather_data
-from preprocessing import scale_data, create_sequences, inverse_scale
 from model_train import (
-    train_linear_regression_gd,
-    predict_linear_regression,
+    cosine_lr_schedule,
+    # cyclical_lr_schedule,
     evaluate_model,
+    # exponential_decay_schedule,
+    # linear_decay_schedule,
+    # polynomial_decay_schedule,
+    predict_linear_regression,
+    run_autoregressive_prediction,
+    # step_decay_schedule,
+    train_linear_regression_gd,
+    warmup_schedule,
+)
+from preprocessing import (
+    add_rolling_features,
+    create_sequences,
+    inverse_scale,
+    scale_data,
 )
 from visualize import (
+    # inspect_residuals,
+    plot_error_hist,  # plot_loss_curve,
+    plot_feature_importance,
+    # plot_interactive_time_series,
+    plot_loss_curve_with_annotations,
     plot_original_vs_scaled,
     plot_predictions,
+    plot_residual_scatter,
     plot_residuals,
-    plot_error_hist,
     plot_scatter,
-    plot_loss_curve,
+    plot_time_series,
     plot_value_distribution,
     visualize_sequence_window,
 )
-from sklearn.model_selection import train_test_split
+
+plt.ioff()
+data_file = "hanoi_weather_data.csv"
+model_file = "theta.npy"
 
 # Step 1: Load data
-data = load_weather_data("hanoi_weather_data.csv")
+data = load_weather_data(data_file)
+
+# Add rolling features
+data = add_rolling_features(data, window=5)
 
 plot_value_distribution(data)
-
+plot_time_series(data, feature_names=["QV2M"], title="QV2M Time Series")
+# plot_interactive_time_series(data, feature_names=["QV2M"])
 # Step 2: Preprocess (scale and sequence)
 scaled_data, scaler = scale_data(data)
 plot_original_vs_scaled(data, scaled_data)
@@ -30,42 +60,127 @@ plot_original_vs_scaled(data, scaled_data)
 seq_length = 50
 X, y = create_sequences(scaled_data, seq_length)
 
+# Flatten X to 2D (n_samples, seq_length * n_features)
+X = X.reshape(X.shape[0], -1)
+
+if y.ndim == 1:
+    y = y.reshape(-1, 1)  # Single target case
+
 visualize_sequence_window(data["QV2M"], seq_length=10, start_idx=1000)
 
 # Step 3: Train/test split
+test_size = 0.2
 date_index = data.index[seq_length:]
 X_train, X_test, y_train, y_test = train_test_split(
-    X, y.reshape(-1), test_size=0.02, shuffle=False
+    X, y, test_size=test_size, shuffle=False
 )
-date_train, date_test = train_test_split(date_index, test_size=0.02, shuffle=False)
+
+date_train, date_test = train_test_split(date_index, test_size=test_size, shuffle=False)
+
+# this is the lr schedulers
+cos_lr = lambda epoch: cosine_lr_schedule(
+    epoch, max_lr=0.005, min_lr=0.0001, total_epochs=500
+)
+# step_lr = lambda epoch: step_decay_schedule(
+#     epoch, initial_lr=0.01, drop=0.5, epochs_drop=100
+# )
+# expo_lr = lambda epoch: exponential_decay_schedule(
+#     epoch, initial_lr=0.01, decay_rate=0.96, decay_steps=100
+# )
+# linear_lr = lambda epoch: linear_decay_schedule(
+#     epoch, initial_lr=0.01, min_lr=0.001, total_epochs=1000
+# )
+# cycle_lr = lambda epoch: cyclical_lr_schedule(
+#     epoch, base_lr=0.001, max_lr=0.01, step_size=200
+# )
+# poly_lr = lambda epoch: polynomial_decay_schedule(
+#     epoch, initial_lr=0.005, min_lr=0.0001, total_epochs=1000, power=2
+# )
+# warmup_lr = lambda epoch: warmup_schedule(
+#     epoch, warmup_epochs=10, initial_lr=0.001, max_lr=0.01
+# )
+
+# Combined learning rate schedule: Warmup + Cosine Annealing
+combined_lr = lambda epoch: (
+    warmup_schedule(epoch, warmup_epochs=10, initial_lr=0.0001, max_lr=0.005)
+    if epoch < 10
+    else cosine_lr_schedule(epoch - 10, max_lr=0.005, min_lr=0.0001, total_epochs=990)
+)
+
+losses = None
+trained = None
+if os.path.isfile(model_file):
+    if True or input("found model, train anyway? ").strip().lower() in "y":
+        print("found model but training")
+        trained = False
+    else:
+        print("found model")
+        trained = True
+else:
+    print("model not found, training")
 
 # Step 4: Train model with gradient descent
-# import math
-# lr_schedule = lambda epoch: max(0.005, 0.01 * (1 + math.cos(epoch / 1000 * math.pi)) / 2)
+if not trained:
+    X_train, X_val, y_train, y_val = train_test_split(
+        X_train, y_train, test_size=0.2, shuffle=False
+    )
 
-theta, losses = train_linear_regression_gd(
-    X_train,
-    y_train,
-    # lr_schedule=lr_schedule,
-    n_epochs=10000,
-)
+    theta, losses, val_losses = train_linear_regression_gd(
+        X_train,
+        y_train,
+        lr_schedule=combined_lr,  # Use the combined learning rate schedule
+        n_epochs=1000,
+        batch_size=64,
+        patience=300,
+        l2_reg=0.001,  # Adjusted regularization strength
+        val_data=(X_val, y_val),  # Pass validation data
+    )
+
+    # Save model
+    np.save(model_file, theta)
+else:
+    # Load model
+    theta = np.load(model_file)
 
 predictions = predict_linear_regression(X_test, theta)
 
 # Step 5: Evaluate
-mse, r2 = evaluate_model(y_test, predictions)
-print(f"MSE: {mse:.4f}, R2: {r2:.4f}")
+mse, r2, mae = evaluate_model(y_test, predictions)
+print(f"MSE: {mse:.4f}, R2: {r2:.4f}, MAE: {mae:.4f}")
 
 # Step 6: Rescale predictions
-y_test_rescaled = inverse_scale(y_test, scaler, data.shape[1])
-predictions_rescaled = inverse_scale(predictions, scaler, data.shape[1])
+y_test_rescaled = inverse_scale(y_test, scaler, y_test.shape[1])
+predictions_rescaled = inverse_scale(predictions, scaler, y_test.shape[1])
+
+# Define feature names
+feature_names = data.columns.tolist()
 
 # Step 7: Visualize
 residuals = y_test - predictions
 errors = abs(residuals)
 
-plot_predictions(date_test, y_test_rescaled, predictions_rescaled)
-plot_residuals(date_test.to_numpy(), residuals)
-plot_error_hist(errors)
+plot_predictions(
+    date_test,
+    y_test_rescaled,
+    predictions_rescaled,
+    feature_names=feature_names,
+    title="Actual vs Predicted",
+)
+plot_residuals(date_test.to_numpy(), residuals, feature_names=feature_names)
+plot_error_hist(errors, feature_names=feature_names)
 plot_scatter(y_test.reshape(-1, 1), predictions.reshape(-1, 1))
-plot_loss_curve(losses)
+plot_residual_scatter(y_test, predictions)
+# == inspect_residuals(predictions, residuals)
+if losses:
+    plot_loss_curve_with_annotations(losses)
+if val_losses:
+    plot_loss_curve_with_annotations(val_losses)
+feature_names = [f"{col}_t-{i}" for col in data.columns for i in range(seq_length)]
+plot_feature_importance(theta, feature_names=feature_names)
+
+# Step 8: Autoregressive Prediction
+autoregressive_predictions = run_autoregressive_prediction(
+    theta, X_test, y_test, scaler, date_test, data.shape[1]
+)
+
+plt.ion()

model_train.py

@@ -1,40 +1,152 @@
+# from itertools import islice
+import time
+
 import numpy as np
-from sklearn.metrics import mean_squared_error, r2_score
+from matplotlib import pyplot as plt
+from sklearn.metrics import (  # type: ignore[import]
+    mean_absolute_error,
+    mean_squared_error,
+    r2_score,
+)
+
+from preprocessing import inverse_scale
 
 
-def train_linear_regression_gd(X_train, y_train, lr_schedule= lambda _e: 0.01, n_epochs=1000):
-    """
-    Train linear regression model using batch gradient descent.
+def train_linear_regression_gd(
+    X_train,
+    y_train,
+    lr_schedule=lambda _e: 0.01,
+    n_epochs=1000,
+    batch_size=32,
+    patience=20,
+    l2_reg=0.0,
+    l1_reg=0.0,
+    val_data=None,  # Add validation data
+):
+    # Ensure X_train is 2D
+    assert X_train.ndim == 2, "X_train must be 2D (n_samples, seq_length * n_features)."
+
+    # Ensure y_train is 2D
+    if y_train.ndim == 1:
+        y_train = y_train.reshape(-1, 1)  # Handle single target case
 
-    Returns:
-        theta: learned weights
-        losses: list of MSE loss per epoch
-    """
     m, n = X_train.shape
-    X_b = np.hstack([np.ones((m, 1)), X_train])
-    theta = np.zeros(n + 1)
+    n_targets = y_train.shape[1]  # Number of target variables
+
+    # Add bias term
+    X_b = np.hstack([np.ones((m, 1)), X_train])  # Add bias term
+
+    # Initialize theta with small random values
+    theta = np.random.randn(n + 1, n_targets) * 0.01
     losses = []
+    val_losses = []  # Track validation losses
+
+    best_loss = float("inf")
+    patience_counter = 0
+
+    X_val, y_val = val_data if val_data else (None, None)
+    if X_val is not None:
+        X_val_b = np.hstack(
+            [np.ones((X_val.shape[0], 1)), X_val]
+        )  # Add bias term to validation data
 
     for epoch in range(n_epochs):
-        lr = lr_schedule(epoch)
-        predictions = X_b @ theta
-        errors = predictions - y_train
-        gradients = (2 / m) * X_b.T @ errors
-        gradients = np.clip(gradients, -1, 1)  # gradient clipping
-        theta -= lr * gradients
-        loss = np.mean(errors ** 2)
-        losses.append(loss)
+        start_time = time.time()  # Start timing the epoch
 
-        if epoch % 1000 == 0:
-            print(f"Epoch {epoch}: Loss = {loss:.6f}, LR: {lr:.6f}")
+        # Shuffle the data at the start of each epoch
+        indices = np.random.permutation(m)
+        X_b_shuffled = X_b[indices]
+        y_shuffled = y_train[indices]
+
+        epoch_loss = 0.0
+        gradient_norm = 0.0
+
+        for i in range(0, m, batch_size):
+            X_batch = X_b_shuffled[i : i + batch_size]
+            y_batch = y_shuffled[i : i + batch_size]
+
+            lr = lr_schedule(epoch)
+            predictions = X_batch @ theta  # Shape: (batch_size, n_targets)
+            errors = predictions - y_batch  # Shape: (batch_size, n_targets)
+
+            # Compute gradients with regularization
+            gradients = (
+                (2 / batch_size) * X_batch.T @ errors
+            )  # Shape: (n_features + 1, n_targets)
+            if l2_reg > 0:
+                gradients[1:] += (
+                    2 * l2_reg * theta[1:]
+                )  # Apply L2 regularization (skip bias)
+            if l1_reg > 0:
+                gradients[1:] += l1_reg * np.sign(
+                    theta[1:]
+                )  # Apply L1 regularization (skip bias)
+
+            # Clip gradients to prevent exploding gradients
+            clip_value = 1e3
+            gradients = np.clip(gradients, -clip_value, clip_value)
+
+            gradient_norm += np.linalg.norm(gradients)  # Track gradient norm
+            theta -= lr * gradients  # Update theta
+
+            batch_loss = np.mean(errors**2)  # Compute batch loss
+            if l2_reg > 0:
+                batch_loss += l2_reg * np.sum(theta[1:] ** 2)  # Add L2 penalty
+            if l1_reg > 0:
+                batch_loss += l1_reg * np.sum(np.abs(theta[1:]))  # Add L1 penalty
+
+            epoch_loss += batch_loss
+
+        epoch_loss /= m // batch_size
+        losses.append(epoch_loss)
+
+        # Compute validation loss
+        if X_val is not None:
+            val_predictions = X_val_b @ theta
+            val_loss = np.mean((val_predictions - y_val) ** 2)
+            val_losses.append(val_loss)
+
+        # Check for early stopping
+        if epoch_loss < best_loss:
+            best_loss = epoch_loss
+            patience_counter = 0
+        else:
+            patience_counter += 1
+
+        if patience_counter >= patience:
+            print(f"Early stopping at epoch {epoch}")
+            break
+
+        # Log statistics
+        epoch_time = time.time() - start_time
+        if epoch % 50 == 0 or epoch == n_epochs - 1:
+            val_loss_str = (
+                f", Val Loss = {val_losses[-1]:.6f}" if X_val is not None else ""
+            )
+            print(
+                f"Epoch {epoch}: Loss = {epoch_loss:.6f}{val_loss_str}, LR: {lr:.6f}, "
+                f"Gradient Norm = {gradient_norm:.4f}, Epoch Time = {epoch_time:.2f}s"
+            )
 
     print(f"Final Loss: {losses[-1]:.6f}")
-    return theta, losses
+    if X_val is not None:
+        print(f"Final Validation Loss: {val_losses[-1]:.6f}")
+        return theta, losses, val_losses
+    return theta, losses, None
 
 
 def predict_linear_regression(X, theta):
-    """Predict using learned linear model."""
-    X_b = np.hstack([np.ones((X.shape[0], 1)), X])
+    """
+    Predict using the linear regression model.
+
+    Args:
+        X: Input features (2D or 3D array).
+        theta: Model weights.
+
+    Returns:
+        predictions: Predicted values.
+    """
+    X_b = np.hstack([np.ones((X.shape[0], 1)), X])  # Add bias term
     return X_b @ theta
 
 
@@ -42,4 +154,197 @@ def evaluate_model(y_true, y_pred):
     """Evaluate using sklearn's built-in MSE and R² functions."""
     mse = mean_squared_error(y_true, y_pred)
     r2 = r2_score(y_true, y_pred)
-    return mse, r2
+    mae = mean_absolute_error(y_true, y_pred)
+    return iter((mse, r2, mae))
+
+
+import math
+
+
+def cosine_lr_schedule(epoch, max_lr=0.01, min_lr=0.001, total_epochs=1000):
+    return (
+        min_lr + (max_lr - min_lr) * (1 + math.cos(math.pi * epoch / total_epochs)) / 2
+    )
+
+
+def step_decay_schedule(epoch, initial_lr=0.01, drop=0.5, epochs_drop=50):
+    """
+    Step decay learning rate scheduler.
+
+    Args:
+        epoch: Current epoch.
+        initial_lr: Initial learning rate.
+        drop: Factor by which to reduce the learning rate.
+        epochs_drop: Number of epochs after which to reduce the learning rate.
+
+    Returns:
+        Updated learning rate.
+    """
+    return initial_lr * (drop ** (epoch // epochs_drop))
+
+
+def exponential_decay_schedule(epoch, initial_lr=0.01, decay_rate=0.1, decay_steps=100):
+    """
+    Exponential decay learning rate scheduler.
+
+    Args:
+        epoch: Current epoch.
+        initial_lr: Initial learning rate.
+        decay_rate: Rate of decay.
+        decay_steps: Number of steps for decay.
+
+    Returns:
+        Updated learning rate.
+    """
+    return initial_lr * (decay_rate ** (epoch / decay_steps))
+
+
+def linear_decay_schedule(epoch, initial_lr=0.01, min_lr=0.001, total_epochs=1000):
+    """
+    Linear decay learning rate scheduler.
+
+    Args:
+        epoch: Current epoch.
+        initial_lr: Initial learning rate.
+        min_lr: Minimum learning rate.
+        total_epochs: Total number of epochs.
+
+    Returns:
+        Updated learning rate.
+    """
+    return max(min_lr, initial_lr - (initial_lr - min_lr) * (epoch / total_epochs))
+
+
+def cyclical_lr_schedule(epoch, base_lr=0.001, max_lr=0.01, step_size=200):
+    """
+    Cyclical learning rate scheduler.
+
+    Args:
+        epoch: Current epoch.
+        base_lr: Minimum learning rate.
+        max_lr: Maximum learning rate.
+        step_size: Number of epochs for half a cycle.
+
+    Returns:
+        Updated learning rate.
+    """
+    cycle = np.floor(1 + epoch / (2 * step_size))
+    x = np.abs(epoch / step_size - 2 * cycle + 1)
+    return base_lr + (max_lr - base_lr) * max(0, (1 - x))
+
+
+def polynomial_decay_schedule(
+    epoch, initial_lr=0.01, min_lr=0.001, total_epochs=1000, power=2
+):
+    """
+    Polynomial decay learning rate scheduler.
+
+    Args:
+        epoch: Current epoch.
+        initial_lr: Initial learning rate.
+        min_lr: Minimum learning rate.
+        total_epochs: Total number of epochs.
+        power: Polynomial power.
+
+    Returns:
+        Updated learning rate.
+    """
+    return min_lr + (initial_lr - min_lr) * ((1 - epoch / total_epochs) ** power)
+
+
+def warmup_schedule(epoch, warmup_epochs=10, initial_lr=0.01, max_lr=0.1):
+    """
+    Warm-up learning rate scheduler.
+
+    Args:
+        epoch: Current epoch.
+        warmup_epochs: Number of warm-up epochs.
+        initial_lr: Initial learning rate.
+        max_lr: Maximum learning rate after warm-up.
+
+    Returns:
+        Updated learning rate.
+    """
+    if epoch < warmup_epochs:
+        return initial_lr + (max_lr - initial_lr) * (epoch / warmup_epochs)
+    return max_lr
+
+
+def autoregressive_prediction(model, initial_sequence, steps, scaler, n_features):
+    """
+    Perform autoregressive predictions using the model.
+
+    Args:
+        model: The trained model (theta for linear regression).
+        initial_sequence: The initial sequence to start predictions (from test data).
+        steps: Number of steps to predict.
+        scaler: Scaler used for inverse scaling.
+        n_features: Number of features in the data.
+
+    Returns:
+        predictions: List of predicted values.
+    """
+    sequence = initial_sequence.copy()  # Start with the initial sequence
+    predictions = []
+
+    for _ in range(steps):
+        # Predict the next value
+        next_prediction = predict_linear_regression(
+            sequence.reshape(1, -1), model
+        )  # Shape: (1, n_targets)
+        predictions.append(next_prediction[0])  # Append the prediction (flattened)
+
+        # Update the sequence: Remove the oldest timestep and append the new prediction
+        sequence = np.roll(sequence, -n_features)  # Shift sequence left by n_features
+        sequence[-n_features:] = (
+            next_prediction  # Add the new prediction (entire target vector)
+        )
+
+    # Rescale predictions back to original scale
+    predictions_rescaled = inverse_scale(np.array(predictions), scaler, n_features)
+    return predictions_rescaled
+
+
+def run_autoregressive_prediction(
+    model, X_test, y_test, scaler, date_test, n_features, steps=None
+):
+    """
+    Run autoregressive prediction and visualize the results.
+
+    Args:
+        model: The trained model (theta for linear regression).
+        X_test: Test features.
+        y_test: Test labels.
+        scaler: Scaler used for inverse scaling.
+        date_test: Dates corresponding to the test set.
+        n_features: Number of features in the data.
+        steps: Number of steps to predict (default is the length of y_test).
+
+    Returns:
+        autoregressive_predictions: Rescaled autoregressive predictions.
+    """
+    if steps is None:
+        steps = len(y_test)  # Default to the length of the test set
+
+    # Use the first test sequence as the initial sequence
+
+    initial_sequence = X_test[0]
+
+    # Perform autoregressive prediction
+    autoregressive_predictions = autoregressive_prediction(
+        model, initial_sequence, steps, scaler, n_features
+    )
+
+    # Visualize the results
+    plt.figure(figsize=(12, 6))
+    plt.plot(date_test, inverse_scale(y_test, scaler, n_features), label="Actual")
+    plt.plot(date_test, autoregressive_predictions, label="Autoregressive Predictions")
+    plt.title("Autoregressive Predictions vs Actual")
+    plt.xlabel("Date")
+    plt.ylabel("Value")
+    plt.legend()
+    plt.grid(True)
+    plt.tight_layout()
+    plt.show()
+
+    return autoregressive_predictions

preprocessing.py

@@ -1,10 +1,26 @@
 import numpy as np
-from sklearn.preprocessing import MinMaxScaler
-from sklearn.model_selection import train_test_split
+import pandas as pd
+from sklearn.preprocessing import MinMaxScaler  # type: ignore[import]
+
+# from sklearn.model_selection import train_test_split
 
 
 def scale_data(data):
-    """Normalize the data using MinMaxScaler."""
+    """Normalize the data using MinMaxScaler.
+
+    Args:
+        data: DataFrame or array with multiple columns.
+
+    Returns:
+        scaled_data: Scaled data as a NumPy array.
+        scaler: Fitted MinMaxScaler object.
+    """
+    if isinstance(data, np.ndarray):
+        if data.ndim == 1:
+            data = data.reshape(-1, 1)  # Ensure 2D for single-column arrays
+    else:
+        data = data.to_numpy()  # Convert DataFrame to NumPy array
+
     scaler = MinMaxScaler()
     scaled_data = scaler.fit_transform(data)
     return scaled_data, scaler
@@ -16,11 +32,44 @@ def inverse_scale(scaled, scaler, n_features):
 
 
 def create_sequences(data, seq_length):
-    """Transform the time series into overlapping sequences."""
-    xs, ys = [], []
+    """
+    Create sequences for time-series data.
+
+    Args:
+        data: Scaled data as a NumPy array (multi-column supported).
+        seq_length: Length of each sequence.
+
+    Returns:
+        X: Input sequences (3D array).
+        y: Target values (2D array for multi-column, 1D for single-column).
+    """
+    X, y = [], []
     for i in range(len(data) - seq_length):
-        x = data[i : i + seq_length].flatten()
-        y = data[i + seq_length]
-        xs.append(x)
-        ys.append(y)
-    return np.array(xs), np.array(ys)
+        X.append(data[i : i + seq_length])
+        y.append(data[i + seq_length])  # Automatically handles multi-column
+    return np.array(X), np.array(y)
+
+
+def add_rolling_features(data, window=5):
+    """
+    Add rolling mean and standard deviation features to the dataset.
+
+    Args:
+        data: DataFrame containing the original features.
+        window: Window size for rolling calculations.
+
+    Returns:
+        DataFrame with additional rolling features.
+    """
+    rolling_data = data.copy()
+    for col in data.columns:
+        rolling_data[f"{col}_rolling_mean_{window}"] = (
+            data[col].rolling(window=window).mean()
+        )
+        rolling_data[f"{col}_rolling_std_{window}"] = (
+            data[col].rolling(window=window).std()
+        )
+
+    # Drop rows with NaN values introduced by rolling calculations
+    rolling_data = rolling_data.dropna()
+    return rolling_data

visualize.py

@@ -2,59 +2,45 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 
-def plot_predictions(date_test, y_test_rescaled, predictions_rescaled):
-    plt.figure(figsize=(12, 6))
-    plt.plot(date_test, y_test_rescaled[:, 0], label="Actual")
-    plt.plot(date_test, predictions_rescaled[:, 0], label="Predicted")
-    plt.title("Actual vs Predicted (First Feature)")
-    plt.xlabel("Date")
-    plt.ylabel("Value")
-    plt.xticks(rotation=45)
-    plt.legend()
-    plt.tight_layout()
-    plt.grid(True)
-    plt.show()
+# data exploration
 
 
-def plot_residuals(date_test, residuals):
-    plt.figure(figsize=(10, 4))
-    plt.plot(date_test, residuals, label="Residuals")
-    plt.hlines(0, 0, len(residuals), colors="red")
-    plt.xlim(date_test.min(), date_test.max())
-    plt.title("Residuals Over Time")
-    plt.xlabel("Date")
-    plt.ylabel("Error")
-    plt.xticks(rotation=45)
-    plt.grid(True)
-    plt.tight_layout()
-    plt.legend()
-    plt.show()
+def plot_value_distribution(data):
+    # import matplotlib.pyplot as plt
 
-
-def plot_error_hist(errors):
-    plt.hist(errors, bins=50)
-    plt.title("Error Distribution")
-    plt.xlabel("Absolute Error")
+    plt.figure(figsize=(8, 4))
+    plt.hist(data.iloc[:, 0], bins=40, color="skyblue", edgecolor="black")
+    plt.title("Distribution of Values")
+    plt.xlabel(data.columns[0])
     plt.ylabel("Frequency")
     plt.grid(True)
+    plt.tight_layout()
     plt.show()
 
 
-def plot_scatter(y_test, predictions):
-    plt.scatter(y_test[:, 0], predictions[:, 0], alpha=0.5)
-    plt.plot([0, 1], [0, 1], "r--")
-    plt.title("Predicted vs Actual (Normalized)")
-    plt.xlabel("Actual")
-    plt.ylabel("Predicted")
-    plt.grid(True)
-    plt.show()
+def plot_time_series(
+    data, feature_names=None, title="Time Series Trends", xlabel="Date", ylabel="Value"
+):
+    """
+    Plot time series data for specified features.
 
+    Args:
+        data: DataFrame containing the time series data.
+        feature_names: List of feature names to plot. If None, all features are plotted.
+        title: Title of the plot.
+        xlabel: Label for the x-axis.
+        ylabel: Label for the y-axis.
+    """
+    if feature_names is None:
+        feature_names = data.columns
 
-def plot_loss_curve(losses):
-    plt.plot(losses)
-    plt.title("Training Loss Curve")
-    plt.xlabel("Epoch")
-    plt.ylabel("MSE Loss")
+    plt.figure(figsize=(12, 6))
+    for feature in feature_names:
+        plt.plot(data.index, data[feature], label=feature)
+    plt.title(title)
+    plt.xlabel(xlabel)
+    plt.ylabel(ylabel)
+    plt.legend()
     plt.grid(True)
     plt.tight_layout()
     plt.show()
@@ -65,7 +51,7 @@ def plot_original_vs_scaled(original_data, scaled_data, feature_names=None):
     Plot original and scaled data, aligned by date.
     """
     if feature_names is None:
-        feature_names = original_data.columns[:2]
+        feature_names = original_data.columns[:3]  # demo
 
     dates = original_data.index  # Grab the datetime index
 
@@ -89,19 +75,6 @@ def plot_original_vs_scaled(original_data, scaled_data, feature_names=None):
     plt.show()
 
 
-def plot_value_distribution(data):
-    # import matplotlib.pyplot as plt
-
-    plt.figure(figsize=(8, 4))
-    plt.hist(data.iloc[:, 0], bins=40, color="skyblue", edgecolor="black")
-    plt.title("Distribution of Values")
-    plt.xlabel(data.columns[0])
-    plt.ylabel("Frequency")
-    plt.grid(True)
-    plt.tight_layout()
-    plt.show()
-
-
 def visualize_sequence_window(series, seq_length=5, start_idx=0):
     """
     Visualize how a sequence window is used to create X (input) and y (target).
@@ -116,18 +89,28 @@ def visualize_sequence_window(series, seq_length=5, start_idx=0):
     else:
         data = series.to_numpy()
 
-    x_seq = data[start_idx:start_idx + seq_length]
+    x_seq = data[start_idx : start_idx + seq_length]
     y_val = data[start_idx + seq_length]
 
     full_seq = np.concatenate([x_seq, [y_val]])
 
     plt.figure(figsize=(8, 3))
-    plt.plot(range(start_idx, start_idx + seq_length + 1), full_seq, marker='o')
+    plt.plot(range(start_idx, start_idx + seq_length + 1), full_seq, marker="o")
     plt.axvline(x=start_idx + seq_length - 0.5, linestyle="--", color="gray")
 
     for i in range(seq_length):
-        plt.text(start_idx + i, x_seq[i], f"x[{i}]", ha="center", va="bottom", fontsize=9)
-    plt.text(start_idx + seq_length, y_val, "y", ha="center", va="bottom", color="red", fontsize=9)
+        plt.text(
+            start_idx + i, x_seq[i], f"x[{i}]", ha="center", va="bottom", fontsize=9
+        )
+    plt.text(
+        start_idx + seq_length,
+        y_val,
+        "y",
+        ha="center",
+        va="bottom",
+        color="red",
+        fontsize=9,
+    )
 
     plt.title(f"Sequence Window (seq_length={seq_length}) → Predict Next")
     plt.xlabel("Time Step")
@@ -135,3 +118,206 @@ def visualize_sequence_window(series, seq_length=5, start_idx=0):
     plt.grid(True)
     plt.tight_layout()
     plt.show()
+
+
+# training analysis
+
+
+def plot_loss_curve_with_annotations(losses, annotate_min=True):
+    """
+    Plot the training loss curve with optional annotations.
+
+    Args:
+        losses: List of loss values for each epoch.
+        annotate_min: Whether to annotate the minimum loss point.
+    """
+    plt.figure(figsize=(10, 6))
+    plt.plot(losses, label="Loss")
+    if annotate_min:
+        min_loss_idx = np.argmin(losses)
+        plt.scatter(
+            min_loss_idx,
+            losses[min_loss_idx],
+            color="red",
+            label=f"Min Loss: {losses[min_loss_idx]:.4f}",
+        )
+    plt.title("Training Loss Curve")
+    plt.xlabel("Epoch")
+    plt.ylabel("MSE Loss")
+    plt.legend()
+    plt.grid(True)
+    plt.tight_layout()
+    plt.show()
+
+
+# post training evaluation
+
+
+def plot_predictions(
+    date_test,
+    y_test_rescaled,
+    predictions_rescaled,
+    feature_names=None,
+    title="Actual vs Predicted",
+    xlabel="Date",
+    ylabel="Value",
+):
+    """
+    Plot actual vs predicted values for each column separately.
+
+    Args:
+        date_test: Array of test dates.
+        y_test_rescaled: Rescaled actual values.
+        predictions_rescaled: Rescaled predicted values.
+        feature_names: List of feature names corresponding to columns.
+        title: Title of the plot.
+        xlabel: Label for the x-axis.
+        ylabel: Label for the y-axis.
+    """
+    n_columns = y_test_rescaled.shape[1]
+    if feature_names is None:
+        feature_names = [f"Column {i}" for i in range(n_columns)]
+
+    for i in range(n_columns):
+        plt.figure(figsize=(12, 6))
+        plt.plot(date_test, y_test_rescaled[:, i], label=f"Actual ({feature_names[i]})")
+        plt.plot(
+            date_test,
+            predictions_rescaled[:, i],
+            label=f"Predicted ({feature_names[i]})",
+        )
+        plt.title(f"{title} ({feature_names[i]})")
+        plt.xlabel(xlabel)
+        plt.ylabel(ylabel)
+        plt.xticks(rotation=45)
+        plt.legend()
+        plt.tight_layout()
+        plt.grid(True)
+        plt.show()
+
+
+def plot_residuals(date_test, residuals, feature_names=None):
+    """
+    Plot residuals over time for each column separately.
+
+    Args:
+        date_test: Array of test dates.
+        residuals: Array of residuals (n_samples, n_columns).
+        feature_names: List of feature names corresponding to columns.
+    """
+    n_columns = residuals.shape[1]
+    if feature_names is None:
+        feature_names = [f"Column {i}" for i in range(n_columns)]
+
+    for i in range(n_columns):
+        plt.figure(figsize=(10, 4))
+        plt.plot(date_test, residuals[:, i], label=f"Residuals ({feature_names[i]})")
+        plt.hlines(0, date_test.min(), date_test.max(), colors="red", linestyles="--")
+        plt.title(f"Residuals Over Time ({feature_names[i]})")
+        plt.xlabel("Date")
+        plt.ylabel("Error")
+        plt.xticks(rotation=45)
+        plt.grid(True)
+        plt.tight_layout()
+        plt.legend()
+        plt.show()
+
+
+def plot_error_hist(errors, feature_names=None):
+    """
+    Plot error histograms for each column separately.
+
+    Args:
+        errors: Array of errors (n_samples, n_columns).
+        feature_names: List of feature names corresponding to columns.
+    """
+    n_columns = errors.shape[1]
+    if feature_names is None:
+        feature_names = [f"Column {i}" for i in range(n_columns)]
+
+    for i in range(n_columns):
+        plt.figure(figsize=(8, 4))
+        plt.hist(errors[:, i], bins=50, alpha=0.7, label=f"{feature_names[i]}")
+        plt.title(f"Error Distribution ({feature_names[i]})")
+        plt.xlabel("Absolute Error")
+        plt.ylabel("Frequency")
+        plt.grid(True)
+        plt.tight_layout()
+        plt.show()
+
+
+def plot_scatter(y_test, predictions):
+    plt.scatter(y_test[:, 0], predictions[:, 0], alpha=0.5)
+    plt.plot([0, 1], [0, 1], "r--")
+    plt.title("Predicted vs Actual (Normalized)")
+    plt.xlabel("Actual")
+    plt.ylabel("Predicted")
+    plt.grid(True)
+    plt.show()
+
+
+def plot_residual_scatter(y_test, predictions):
+    residuals = y_test - predictions
+    plt.figure(figsize=(8, 6))
+    plt.scatter(predictions, residuals, alpha=0.5)
+    plt.axhline(0, color="red", linestyle="--")
+    plt.title("Residuals vs Predictions")
+    plt.xlabel("Predicted Values")
+    plt.ylabel("Residuals")
+    plt.grid(True)
+    plt.tight_layout()
+    plt.show()
+
+
+def plot_feature_importance(theta, feature_names):
+    """
+    Plot feature importance based on learned weights.
+
+    Args:
+        theta: Array of weights (including bias) with shape (n_features + 1, n_targets).
+        feature_names: List of feature names.
+    """
+    if feature_names is None:
+        feature_names = [f"Feature {i}" for i in range(theta.shape[0] - 1)]
+
+    if theta.shape[0] != len(feature_names) + 1:
+        raise ValueError(
+            f"Mismatch: theta has {theta.shape[0]} rows (including bias), "
+            f"but feature_names has {len(feature_names)} elements."
+        )
+
+    # Average weights across all targets
+    avg_theta = np.mean(theta, axis=1)
+
+    plt.figure(figsize=(10, 6))
+    plt.bar(range(len(avg_theta)), avg_theta, tick_label=["Bias"] + feature_names)
+    plt.xticks(rotation=45, ha="right")  # Rotate labels 45 degrees
+    plt.title("Feature Importance (Averaged Across Targets)")
+    plt.xlabel("Features")
+    plt.ylabel("Weight")
+    plt.grid(True)
+    plt.tight_layout()
+    plt.show()
+
+
+# import plotly.express as px
+
+
+# def plot_interactive_time_series(data, feature_names=None):
+#     if feature_names is None:
+#         feature_names = data.columns
+
+#     fig = px.line(data, x=data.index, y=feature_names, title="Interactive Time Series")
+#     fig.show()
+
+
+def inspect_residuals(predictions, residuals):
+    plt.figure(figsize=(10, 6))
+    plt.scatter(predictions, residuals, alpha=0.5)
+    plt.axhline(0, color="red", linestyle="--")
+    plt.title("Residuals vs Predictions")
+    plt.xlabel("Predicted Values")
+    plt.ylabel("Residuals")
+    plt.grid(True)
+    plt.tight_layout()
+    plt.show()