this thing runs somehow

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 = {

main.py

@@ -52,14 +52,21 @@ 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.025
+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=test_size, shuffle=False
+    X, y, test_size=test_size, shuffle=False
 )
+
 date_train, date_test = train_test_split(date_index, test_size=test_size, shuffle=False)
 
 # this is the lr schedulers
@@ -82,9 +89,10 @@ cycle_lr = lambda epoch: cyclical_lr_schedule(
 # warmup_lr = lambda epoch: warmup_schedule(
 #     epoch, warmup_epochs=10, initial_lr=0.001, max_lr=0.01
 # )
+losses = None
 trained = None
 if os.path.isfile(model_file):
-    if input("found model, train anyway? ").strip().lower() in "y":
+    if True or input("found model, train anyway? ").strip().lower() in "y":
         print("found model but training")
         trained = False
     else:
@@ -117,15 +125,15 @@ 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])
 
 # Step 7: Visualize
 residuals = y_test - predictions
 errors = abs(residuals)
 
 plot_predictions(
-    date_test, y_test_rescaled, predictions_rescaled, title="Actual vs Predicted (QV2M)"
+    date_test, y_test_rescaled, predictions_rescaled, title="Actual vs Predicted"
 )
 plot_residuals(date_test.to_numpy(), residuals)
 plot_error_hist(errors)
@@ -134,8 +142,8 @@ plot_residual_scatter(y_test, predictions)
 # == inspect_residuals(predictions, residuals)
 if losses:
     plot_loss_curve_with_annotations(losses)
-plot_feature_importance(theta, feature_names=[f"QV2M_t-{i}" for i in range(seq_length)])
-
+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(

model_train.py

@@ -3,8 +3,11 @@ import time
 
 import numpy as np
 from matplotlib import pyplot as plt
-from sklearn.metrics import (mean_absolute_error, mean_squared_error,  # type: ignore[import]
-                             r2_score)
+from sklearn.metrics import (  # type: ignore[import]
+    mean_absolute_error,
+    mean_squared_error,
+    r2_score,
+)
 
 from preprocessing import inverse_scale
 
@@ -19,9 +22,21 @@ def train_linear_regression_gd(
     l2_reg=0.0,
     l1_reg=0.0,
 ):
+    # 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
+
     m, n = X_train.shape
+    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
-    theta = np.zeros(n + 1)
+
+    # Initialize theta
+    theta = np.zeros((n + 1, n_targets))  # (n_features + 1, n_targets)
     losses = []
 
     best_loss = float("inf")
@@ -43,33 +58,36 @@ def train_linear_regression_gd(
             y_batch = y_shuffled[i : i + batch_size]
 
             lr = lr_schedule(epoch)
-            predictions = X_batch @ theta
-            errors = predictions - y_batch
+            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
+            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:]
+                gradients[1:] += (
+                    2 * l2_reg * theta[1:]
+                )  # Apply L2 regularization (skip bias)
             if l1_reg > 0:
-                gradients[1:] += l1_reg * np.sign(theta[1:])
+                gradients[1:] += l1_reg * np.sign(
+                    theta[1:]
+                )  # Apply L1 regularization (skip bias)
 
             gradient_norm += np.linalg.norm(gradients)  # Track gradient norm
-            theta -= lr * gradients
+            theta -= lr * gradients  # Update theta
 
-            batch_loss = np.mean(errors**2)
+            batch_loss = np.mean(errors**2)  # Compute batch loss
             if l2_reg > 0:
-                batch_loss += l2_reg * np.sum(theta[1:] ** 2)
+                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:]))
+                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 weight norm
-        weight_norm = np.linalg.norm(theta)
-
         # Check for early stopping
         if epoch_loss < best_loss:
             best_loss = epoch_loss
@@ -86,8 +104,7 @@ def train_linear_regression_gd(
         if epoch % 50 == 0 or epoch == n_epochs - 1:
             print(
                 f"Epoch {epoch}: Loss = {epoch_loss:.6f}, LR: {lr:.6f}, "
-                f"Gradient Norm = {gradient_norm:.4f}, Weight Norm = {weight_norm:.4f}, "
-                f"Epoch Time = {epoch_time:.2f}s"
+                f"Gradient Norm = {gradient_norm:.4f}, Epoch Time = {epoch_time:.2f}s"
             )
 
     print(f"Final Loss: {losses[-1]:.6f}")
@@ -95,8 +112,17 @@ def train_linear_regression_gd(
 
 
 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
 
 
@@ -239,12 +265,16 @@ def autoregressive_prediction(model, initial_sequence, steps, scaler, n_features
 
     for _ in range(steps):
         # Predict the next value
-        next_prediction = predict_linear_regression(sequence.reshape(1, -1), model)
-        predictions.append(next_prediction[0])
+        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 value and append the new prediction
-        sequence = np.roll(sequence, -1)  # Shift sequence left
-        sequence[-1] = next_prediction  # Add the new prediction
+        # 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)

preprocessing.py

@@ -3,9 +3,24 @@ 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.
+    """
+    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
+    else:
+        data = data.to_numpy()  # Convert DataFrame to NumPy array
+
     scaler = MinMaxScaler()
     scaled_data = scaler.fit_transform(data)
     return scaled_data, scaler
@@ -17,11 +32,19 @@ 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)

visualize.py

@@ -51,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
 
@@ -162,7 +162,7 @@ def plot_predictions(
     ylabel="Value",
 ):
     """
-    Plot actual vs predicted values.
+    Plot actual vs predicted values for each column separately.
 
     Args:
         date_test: Array of test dates.
@@ -172,41 +172,61 @@ def plot_predictions(
         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(title)
-    plt.xlabel(xlabel)
-    plt.ylabel(ylabel)
-    plt.xticks(rotation=45)
-    plt.legend()
-    plt.tight_layout()
-    plt.grid(True)
-    plt.show()
+    n_columns = y_test_rescaled.shape[1]
+    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, predictions_rescaled[:, i], label=f"Predicted (Column {i})")
+        plt.title(f"{title} (Column {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):
-    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()
+    """
+    Plot residuals over time for each column separately.
+
+    Args:
+        date_test: Array of test dates.
+        residuals: Array of residuals (n_samples, n_columns).
+    """
+    n_columns = residuals.shape[1]
+    for i in range(n_columns):
+        plt.figure(figsize=(10, 4))
+        plt.plot(date_test, residuals[:, i], label=f"Residuals (Column {i})")
+        plt.hlines(0, date_test.min(), date_test.max(), colors="red", linestyles="--")
+        plt.title(f"Residuals Over Time (Column {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):
-    plt.hist(errors, bins=50)
-    plt.title("Error Distribution")
-    plt.xlabel("Absolute Error")
-    plt.ylabel("Frequency")
-    plt.grid(True)
-    plt.show()
+    """
+    Plot error histograms for each column separately.
+
+    Args:
+        errors: Array of errors (n_samples, n_columns).
+    """
+    n_columns = errors.shape[1]
+    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.title(f"Error Distribution (Column {i})")
+        plt.xlabel("Absolute Error")
+        plt.ylabel("Frequency")
+        plt.grid(True)
+        plt.tight_layout()
+        plt.show()
 
 
 def plot_scatter(y_test, predictions):
@@ -234,30 +254,36 @@ def plot_residual_scatter(y_test, predictions):
 
 def plot_feature_importance(theta, feature_names):
     """
-    Plot feature importance based on learned weights.
+    Plot feature importance based on learned weights for each target variable.
 
     Args:
-        theta: Array of weights (including bias).
+        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(len(theta) - 1)]
+        feature_names = [f"Feature {i}" for i in range(theta.shape[0] - 1)]
 
-    if len(theta) != len(feature_names) + 1:
+    if theta.shape[0] != len(feature_names) + 1:
         raise ValueError(
-            f"Mismatch: theta has {len(theta)} elements (including bias), "
+            f"Mismatch: theta has {theta.shape[0]} rows (including bias), "
             f"but feature_names has {len(feature_names)} elements."
         )
 
-    plt.figure(figsize=(10, 6))
-    plt.bar(range(len(theta)), theta, tick_label=["Bias"] + feature_names)
-    plt.xticks(rotation=45, ha="right")  # Rotate labels 45 degrees
-    plt.title("Feature Importance (Learned Weights)")
-    plt.xlabel("Features")
-    plt.ylabel("Weight")
-    plt.grid(True)
-    plt.tight_layout()
-    plt.show()
+    n_targets = theta.shape[1]  # Number of target variables
+    for target_idx in range(n_targets):
+        plt.figure(figsize=(10, 6))
+        plt.bar(
+            range(theta.shape[0]),
+            theta[:, target_idx],
+            tick_label=["Bias"] + feature_names,
+        )
+        plt.xticks(rotation=45, ha="right")  # Rotate labels 45 degrees
+        plt.title(f"Feature Importance (Target {target_idx})")
+        plt.xlabel("Features")
+        plt.ylabel("Weight")
+        plt.grid(True)
+        plt.tight_layout()
+        plt.show()
 
 
 # import plotly.express as px
@@ -280,4 +306,4 @@ def inspect_residuals(predictions, residuals):
     plt.ylabel("Residuals")
     plt.grid(True)
     plt.tight_layout()
-    plt.show()
+    plt.show()