real

main.py

@@ -6,8 +6,8 @@ from sklearn.model_selection import train_test_split  # type: ignore[import]
 
 from data_loader import load_weather_data
 from model_train import (
-    # cosine_lr_schedule,
+    cosine_lr_schedule,
-    cyclical_lr_schedule,
+    # cyclical_lr_schedule,
     evaluate_model,
     # exponential_decay_schedule,
     # linear_decay_schedule,
@@ -16,9 +16,14 @@ from model_train import (
     run_autoregressive_prediction,
     # step_decay_schedule,
     train_linear_regression_gd,
-    # warmup_schedule,
+    warmup_schedule,
+)
+from preprocessing import (
+    add_rolling_features,
+    create_sequences,
+    inverse_scale,
+    scale_data,
 )
-from preprocessing import create_sequences, inverse_scale, scale_data
 from visualize import (
     # inspect_residuals,
     plot_error_hist,  # plot_loss_curve,
@@ -42,6 +47,9 @@ model_file = "theta.npy"
 # Step 1: Load data
 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"])
@@ -70,7 +78,9 @@ X_train, X_test, y_train, y_test = train_test_split(
 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.01, min_lr=0.001)
+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
 # )
@@ -80,15 +90,23 @@ date_train, date_test = train_test_split(date_index, test_size=test_size, shuffl
 # 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(
+# cycle_lr = lambda epoch: cyclical_lr_schedule(
-    epoch, base_lr=0.001, max_lr=0.01, step_size=200
+#     epoch, base_lr=0.001, max_lr=0.01, step_size=200
-)
+# )
 # poly_lr = lambda epoch: polynomial_decay_schedule(
-#     epoch, initial_lr=0.01, min_lr=0.001, total_epochs=1000, power=2
+#     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):
@@ -103,13 +121,19 @@ else:
 
 # Step 4: Train model with gradient descent
 if not trained:
-    theta, losses = train_linear_regression_gd(
+    X_train, X_val, y_train, y_val = train_test_split(
+        X_train, y_train, test_size=0.2, shuffle=False
+    )
+
+    theta, losses, val_losses = train_linear_regression_gd(
         X_train,
         y_train,
-        lr_schedule=cycle_lr,
+        lr_schedule=combined_lr,  # Use the combined learning rate schedule
         n_epochs=1000,
-        batch_size=32,
+        batch_size=64,
-        patience=50,
+        patience=300,
+        l2_reg=0.001,  # Adjusted regularization strength
+        val_data=(X_val, y_val),  # Pass validation data
     )
 
     # Save model
@@ -128,20 +152,29 @@ print(f"MSE: {mse:.4f}, R2: {r2:.4f}, MAE: {mae:.4f}")
 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, title="Actual vs Predicted"
+    date_test,
+    y_test_rescaled,
+    predictions_rescaled,
+    feature_names=feature_names,
+    title="Actual vs Predicted",
 )
-plot_residuals(date_test.to_numpy(), residuals)
+plot_residuals(date_test.to_numpy(), residuals, feature_names=feature_names)
-plot_error_hist(errors)
+plot_error_hist(errors, feature_names=feature_names)
 plot_scatter(y_test.reshape(-1, 1), predictions.reshape(-1, 1))
 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)
 

model_train.py

@@ -21,6 +21,7 @@ def train_linear_regression_gd(
     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)."
@@ -35,13 +36,20 @@ def train_linear_regression_gd(
     # Add bias term
     X_b = np.hstack([np.ones((m, 1)), X_train])  # Add bias term
 
-    # Initialize theta
+    # Initialize theta with small random values
-    theta = np.zeros((n + 1, n_targets))  # (n_features + 1, n_targets)
+    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):
         start_time = time.time()  # Start timing the epoch
 
@@ -74,6 +82,10 @@ def train_linear_regression_gd(
                     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
 
@@ -88,6 +100,12 @@ def train_linear_regression_gd(
         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
@@ -102,13 +120,19 @@ def train_linear_regression_gd(
         # 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}, LR: {lr:.6f}, "
+                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):
@@ -303,6 +327,7 @@ def run_autoregressive_prediction(
         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

preprocessing.py

@@ -1,8 +1,10 @@
 import numpy as np
+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.
 
@@ -13,8 +15,6 @@ def scale_data(data):
         scaled_data: Scaled data as a NumPy array.
         scaler: Fitted MinMaxScaler object.
     """
-    from sklearn.preprocessing import MinMaxScaler
-
     if isinstance(data, np.ndarray):
         if data.ndim == 1:
             data = data.reshape(-1, 1)  # Ensure 2D for single-column arrays
@@ -48,3 +48,28 @@ def create_sequences(data, seq_length):
         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

@@ -157,6 +157,7 @@ def plot_predictions(
     date_test,
     y_test_rescaled,
     predictions_rescaled,
+    feature_names=None,
     title="Actual vs Predicted",
     xlabel="Date",
     ylabel="Value",
@@ -168,16 +169,24 @@ def plot_predictions(
         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 (Column {i})")
+        plt.plot(date_test, y_test_rescaled[:, i], label=f"Actual ({feature_names[i]})")
-        plt.plot(date_test, predictions_rescaled[:, i], label=f"Predicted (Column {i})")
+        plt.plot(
-        plt.title(f"{title} (Column {i})")
+            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)
@@ -187,20 +196,24 @@ def plot_predictions(
         plt.show()
 
 
-def plot_residuals(date_test, residuals):
+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 (Column {i})")
+        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 (Column {i})")
+        plt.title(f"Residuals Over Time ({feature_names[i]})")
         plt.xlabel("Date")
         plt.ylabel("Error")
         plt.xticks(rotation=45)
@@ -210,18 +223,22 @@ def plot_residuals(date_test, residuals):
         plt.show()
 
 
-def plot_error_hist(errors):
+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"Column {i}")
+        plt.hist(errors[:, i], bins=50, alpha=0.7, label=f"{feature_names[i]}")
-        plt.title(f"Error Distribution (Column {i})")
+        plt.title(f"Error Distribution ({feature_names[i]})")
         plt.xlabel("Absolute Error")
         plt.ylabel("Frequency")
         plt.grid(True)
@@ -254,7 +271,7 @@ def plot_residual_scatter(y_test, predictions):
 
 def plot_feature_importance(theta, feature_names):
     """
-    Plot feature importance based on learned weights for each target variable.
+    Plot feature importance based on learned weights.
 
     Args:
         theta: Array of weights (including bias) with shape (n_features + 1, n_targets).
@@ -269,21 +286,18 @@ def plot_feature_importance(theta, feature_names):
             f"but feature_names has {len(feature_names)} elements."
         )
 
-    n_targets = theta.shape[1]  # Number of target variables
+    # Average weights across all targets
-    for target_idx in range(n_targets):
+    avg_theta = np.mean(theta, axis=1)
-        plt.figure(figsize=(10, 6))
+
-        plt.bar(
+    plt.figure(figsize=(10, 6))
-            range(theta.shape[0]),
+    plt.bar(range(len(avg_theta)), avg_theta, tick_label=["Bias"] + feature_names)
-            theta[:, target_idx],
+    plt.xticks(rotation=45, ha="right")  # Rotate labels 45 degrees
-            tick_label=["Bias"] + feature_names,
+    plt.title("Feature Importance (Averaged Across Targets)")
-        )
+    plt.xlabel("Features")
-        plt.xticks(rotation=45, ha="right")  # Rotate labels 45 degrees
+    plt.ylabel("Weight")
-        plt.title(f"Feature Importance (Target {target_idx})")
+    plt.grid(True)
-        plt.xlabel("Features")
+    plt.tight_layout()
-        plt.ylabel("Weight")
+    plt.show()
-        plt.grid(True)
-        plt.tight_layout()
-        plt.show()
 
 
 # import plotly.express as px