WAY BETTER

main.py

@@ -1,28 +1,35 @@
 # import pandas as pd
+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,
+    cosine_lr_schedule,
-    predict_linear_regression,
     evaluate_model,
+    predict_linear_regression,
+    train_linear_regression_gd,
 )
+from preprocessing import create_sequences, inverse_scale, scale_data
 from visualize import (
+    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
 
 # Step 1: Load data
 data = load_weather_data("hanoi_weather_data.csv")
 
 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)
@@ -33,20 +40,20 @@ X, y = create_sequences(scaled_data, seq_length)
 visualize_sequence_window(data["QV2M"], seq_length=10, start_idx=1000)
 
 # Step 3: Train/test split
+test_size = 0.025
 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.reshape(-1), 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)
 
 # 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)
 
 theta, losses = train_linear_regression_gd(
     X_train,
     y_train,
-    # lr_schedule=lr_schedule,
+    lr_schedule=lambda epoch: cosine_lr_schedule(epoch, max_lr=0.01, min_lr=0.001),
     # n_epochs=1000,
 )
 
@@ -64,8 +71,12 @@ predictions_rescaled = inverse_scale(predictions, scaler, data.shape[1])
 residuals = y_test - predictions
 errors = abs(residuals)
 
-plot_predictions(date_test, y_test_rescaled, predictions_rescaled)
+plot_predictions(
+    date_test, y_test_rescaled, predictions_rescaled, title="Actual vs Predicted (QV2M)"
+)
 plot_residuals(date_test.to_numpy(), residuals)
 plot_error_hist(errors)
 plot_scatter(y_test.reshape(-1, 1), predictions.reshape(-1, 1))
-plot_loss_curve(losses)
+plot_residual_scatter(y_test, predictions)
+plot_loss_curve_with_annotations(losses)
+plot_feature_importance(theta, feature_names=[f"QV2M_t-{i}" for i in range(seq_length)])

model_train.py

@@ -1,9 +1,16 @@
-from itertools import islice
+# from itertools import islice
 import numpy as np
-from sklearn.metrics import mean_squared_error, r2_score
+from sklearn.metrics import mean_squared_error, r2_score # type: ignore[import]
 
 
-def train_linear_regression_gd(X_train, y_train, lr_schedule=lambda _e: 0.01, n_epochs=1000, batch_size=32):
+def train_linear_regression_gd(
+    X_train,
+    y_train,
+    lr_schedule=lambda _e: 0.01,
+    n_epochs=1000,
+    batch_size=32,
+    patience=20,
+):
     """
     Train linear regression model using batch gradient descent.
 
@@ -23,6 +30,9 @@ def train_linear_regression_gd(X_train, y_train, lr_schedule=lambda _e: 0.01, n_
     theta = np.zeros(n + 1)
     losses = []
 
+    best_loss = float("inf")
+    patience_counter = 0
+
     for epoch in range(n_epochs):
         # Shuffle the data at the start of each epoch
         indices = np.random.permutation(m)
@@ -31,8 +41,8 @@ def train_linear_regression_gd(X_train, y_train, lr_schedule=lambda _e: 0.01, n_
 
         epoch_loss = 0.0
         for i in range(0, m, batch_size):
-            X_batch = X_b_shuffled[i:i + batch_size]
+            X_batch = X_b_shuffled[i : i + batch_size]
-            y_batch = y_shuffled[i:i + batch_size]
+            y_batch = y_shuffled[i : i + batch_size]
 
             lr = lr_schedule(epoch)
             predictions = X_batch @ theta
@@ -41,10 +51,21 @@ def train_linear_regression_gd(X_train, y_train, lr_schedule=lambda _e: 0.01, n_
             gradients = np.clip(gradients, -1, 1)  # Gradient clipping
             theta -= lr * gradients
 
-            batch_loss = np.mean(errors ** 2)
+            batch_loss = np.mean(errors**2)
             epoch_loss += batch_loss
 
-        losses.append(epoch_loss / (m // batch_size))
+        epoch_loss /= m // batch_size
+        losses.append(epoch_loss)
+
+        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
 
         if epoch % 100 == 0:
             print(f"Epoch {epoch}: Loss = {losses[-1]:.6f}, LR: {lr:.6f}")
@@ -63,4 +84,13 @@ 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
+    return mse, r2
+
+
+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
+    )

preprocessing.py

@@ -1,6 +1,6 @@
 import numpy as np
-from sklearn.preprocessing import MinMaxScaler
+from sklearn.preprocessing import MinMaxScaler # type: ignore[import]
-from sklearn.model_selection import train_test_split
+# from sklearn.model_selection import train_test_split
 
 
 def scale_data(data):

visualize.py

@@ -2,13 +2,182 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 
-def plot_predictions(date_test, y_test_rescaled, predictions_rescaled):
+# data exploration
+
+
+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 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
+
+    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()
+
+
+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]
+
+    dates = original_data.index  # Grab the datetime index
+
+    fig, axs = plt.subplots(2, 1, figsize=(12, 6), sharex=True)
+
+    for col in feature_names:
+        axs[0].plot(dates, original_data[col], label=f"Original {col}")
+    axs[0].set_title("Original Data")
+    axs[0].legend()
+    axs[0].grid(True)
+
+    for i, col in enumerate(feature_names):
+        col_index = original_data.columns.get_loc(col)
+        axs[1].plot(dates, scaled_data[:, col_index], label=f"Scaled {col}")
+    axs[1].set_title("Scaled Data (MinMaxScaler)")
+    axs[1].legend()
+    axs[1].grid(True)
+
+    plt.xlabel("Date")
+    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).
+
+    Parameters:
+        series (pd.Series or np.ndarray): 1D time series data
+        seq_length (int): Length of input window
+        start_idx (int): Where to start slicing from for demo
+    """
+    if isinstance(series, np.ndarray):
+        data = series
+    else:
+        data = series.to_numpy()
+
+    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.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.title(f"Sequence Window (seq_length={seq_length}) → Predict Next")
+    plt.xlabel("Time Step")
+    plt.ylabel("Value")
+    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,
+    title="Actual vs Predicted",
+    xlabel="Date",
+    ylabel="Value",
+):
+    """
+    Plot actual vs predicted values.
+
+    Args:
+        date_test: Array of test dates.
+        y_test_rescaled: Rescaled actual values.
+        predictions_rescaled: Rescaled predicted values.
+        title: Title of the plot.
+        xlabel: Label for the x-axis.
+        ylabel: Label for the y-axis.
+    """
     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.title(title)
-    plt.xlabel("Date")
+    plt.xlabel(xlabel)
-    plt.ylabel("Value")
+    plt.ylabel(ylabel)
     plt.xticks(rotation=45)
     plt.legend()
     plt.tight_layout()
@@ -50,88 +219,53 @@ def plot_scatter(y_test, predictions):
     plt.show()
 
 
-def plot_loss_curve(losses):
+def plot_residual_scatter(y_test, predictions):
-    plt.plot(losses)
+    residuals = y_test - predictions
-    plt.title("Training Loss Curve")
+    plt.figure(figsize=(8, 6))
-    plt.xlabel("Epoch")
+    plt.scatter(predictions, residuals, alpha=0.5)
-    plt.ylabel("MSE Loss")
+    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_original_vs_scaled(original_data, scaled_data, feature_names=None):
+def plot_feature_importance(theta, feature_names):
     """
-    Plot original and scaled data, aligned by date.
+    Plot feature importance based on learned weights.
+
+    Args:
+        theta: Array of weights (including bias).
+        feature_names: List of feature names.
     """
     if feature_names is None:
-        feature_names = original_data.columns[:2]
+        feature_names = [f"Feature {i}" for i in range(len(theta) - 1)]
 
-    dates = original_data.index  # Grab the datetime index
+    if len(theta) != len(feature_names) + 1:
+        raise ValueError(
+            f"Mismatch: theta has {len(theta)} elements (including bias), "
+            f"but feature_names has {len(feature_names)} elements."
+        )
 
-    fig, axs = plt.subplots(2, 1, figsize=(12, 6), sharex=True)
+    plt.figure(figsize=(10, 6))
-
+    plt.bar(range(len(theta)), theta, tick_label=["Bias"] + feature_names)
-    for col in feature_names:
+    plt.xticks(rotation=45, ha="right")  # Rotate labels 45 degrees
-        axs[0].plot(dates, original_data[col], label=f"Original {col}")
+    plt.title("Feature Importance (Learned Weights)")
-    axs[0].set_title("Original Data")
+    plt.xlabel("Features")
-    axs[0].legend()
+    plt.ylabel("Weight")
-    axs[0].grid(True)
-
-    for i, col in enumerate(feature_names):
-        col_index = original_data.columns.get_loc(col)
-        axs[1].plot(dates, scaled_data[:, col_index], label=f"Scaled {col}")
-    axs[1].set_title("Scaled Data (MinMaxScaler)")
-    axs[1].legend()
-    axs[1].grid(True)
-
-    plt.xlabel("Date")
-    plt.tight_layout()
-    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):
+# import plotly.express as px
-    """
-    Visualize how a sequence window is used to create X (input) and y (target).
 
-    Parameters:
-        series (pd.Series or np.ndarray): 1D time series data
-        seq_length (int): Length of input window
-        start_idx (int): Where to start slicing from for demo
-    """
-    if isinstance(series, np.ndarray):
-        data = series
-    else:
-        data = series.to_numpy()
 
-    x_seq = data[start_idx:start_idx + seq_length]
+# def plot_interactive_time_series(data, feature_names=None):
-    y_val = data[start_idx + seq_length]
+#     if feature_names is None:
+#         feature_names = data.columns
 
-    full_seq = np.concatenate([x_seq, [y_val]])
+#     fig = px.line(data, x=data.index, y=feature_names, title="Interactive Time Series")
-
+#     fig.show()
-    plt.figure(figsize=(8, 3))
-    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.title(f"Sequence Window (seq_length={seq_length}) → Predict Next")
-    plt.xlabel("Time Step")
-    plt.ylabel("Value")
-    plt.grid(True)
-    plt.tight_layout()
-    plt.show()