real

data_loader.py

@@ -0,0 +1,43 @@
+import pandas as pd
+import io
+from typing import Union
+
+
+def cut_headers(stream: Union[io.StringIO, io.FileIO, io.TextIOWrapper]):
+    """Remove custom header metadata from the CSV file."""
+    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
+
+
+def load_weather_data(file_path: str) -> pd.DataFrame:
+    """Load weather data and preprocess the datetime index."""
+    with open(file_path, "r", encoding="utf-8") as file:
+        _, csv_file = cut_headers(file)
+        data = pd.read_csv(csv_file, na_values=-999)
+
+    # 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"])}
+        )
+        data.drop(columns=["YEAR", "MO", "DY"], inplace=True)
+        data.set_index("date", inplace=True)
+
+    data.dropna(inplace=True)
+    data = data.asfreq("D")  # Assumes regular daily data
+    return data

dataset_update.py

@@ -1,11 +1,12 @@
 import urllib.parse
 import urllib.request
 from datetime import datetime
-from global_config import DATA_FILE_NAME
+
+data_file = "hanoi_weather_data.csv"
 
 # Define parameters
 latitude, longitude = 21.04815477090267, 105.80123159165096  # lmao
-start_date = "2000-01-01"
+start_date = "1985-01-01"
 end_date = datetime.now().strftime("%Y-%m-%d")  # Today's date
 
 # Define desired parameters
@@ -18,9 +19,7 @@ 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 = [
-
     "QV2M",
-
 ]
 
 params = {
@@ -42,13 +41,14 @@ base_url = "https://power.larc.nasa.gov/api/temporal/daily/point"
 url = f"{base_url}?{query_string}"
 
 # Make the API request using urllib
+if __name__ == "__main__":
     try:
         with urllib.request.urlopen(url) as response:
             if response.status == 200:
-            csv_data = response.read().decode("utf-8")
+                csv_data = response.read()
-            with open(DATA_FILE_NAME, "w", encoding="utf-8") as file:
+                with open(data_file, "wb") as file:
                     file.write(csv_data)
-            print("Data successfully downloaded and saved to 'hanoi_weather_data.csv'.")
+                print(f"Data successfully downloaded and saved to {data_file}.")
             else:
                 print(f"Failed to retrieve data: HTTP Status Code {response.status}")
     except urllib.error.URLError as e:

main.py

@@ -0,0 +1,53 @@
+import pandas as pd
+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,
+    evaluate_model,
+)
+from visualize import (
+    plot_predictions,
+    plot_residuals,
+    plot_error_hist,
+    plot_scatter,
+    plot_loss_curve,
+)
+from sklearn.model_selection import train_test_split
+
+# Step 1: Load data
+data = load_weather_data("hanoi_weather_data.csv")
+
+# Step 2: Preprocess (scale and sequence)
+scaled_data, scaler = scale_data(data)
+seq_length = 50
+X, y = create_sequences(scaled_data, seq_length)
+
+# Step 3: Train/test split
+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
+)
+date_train, date_test = train_test_split(date_index, test_size=0.02, shuffle=False)
+
+# Step 4: Train model with gradient descent
+theta, losses = train_linear_regression_gd(X_train, y_train, lr=0.005, n_epochs=300)
+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}")
+
+# Step 6: Rescale predictions
+y_test_rescaled = inverse_scale(y_test, scaler, data.shape[1])
+predictions_rescaled = inverse_scale(predictions, scaler, data.shape[1])
+
+# Step 7: Visualize
+residuals = y_test - predictions
+errors = abs(residuals)
+
+plot_predictions(date_test, y_test_rescaled, predictions_rescaled)
+plot_residuals(date_test, residuals)
+plot_error_hist(errors)
+plot_scatter(y_test.reshape(-1, 1), predictions.reshape(-1, 1))
+plot_loss_curve(losses)

model_train.py

@@ -0,0 +1,45 @@
+import numpy as np
+
+
+def train_linear_regression_gd(X_train, y_train, lr=0.005, n_epochs=300):
+    """
+    Train linear regression model using batch gradient descent.
+
+    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)
+    losses = []
+
+    for epoch in range(n_epochs):
+        predictions = X_b @ theta
+        errors = predictions - y_train
+        gradients = (2 / m) * X_b.T @ errors
+        gradients = np.clip(gradients, -1, 1)  # gradient clipping to prevent divergence
+        theta -= lr * gradients
+        loss = np.mean(errors ** 2)
+        losses.append(loss)
+
+        # Debug print every 10 epochs
+        if epoch % 10 == 0:
+            print(f"Epoch {epoch}: Loss = {loss:.6f}")
+
+    return theta, losses
+
+
+def predict_linear_regression(X, theta):
+    """Predict using learned linear model."""
+    X_b = np.hstack([np.ones((X.shape[0], 1)), X])
+    return X_b @ theta
+
+
+def evaluate_model(y_true, y_pred):
+    """Evaluate MSE and R^2 manually."""
+    mse = np.mean((y_true - y_pred) ** 2)
+    ss_total = np.sum((y_true - np.mean(y_true)) ** 2)
+    ss_res = np.sum((y_true - y_pred) ** 2)
+    r2 = 1 - (ss_res / ss_total)
+    return mse, r2

preprocessing.py

@@ -0,0 +1,26 @@
+import numpy as np
+from sklearn.preprocessing import MinMaxScaler
+from sklearn.model_selection import train_test_split
+
+
+def scale_data(data):
+    """Normalize the data using MinMaxScaler."""
+    scaler = MinMaxScaler()
+    scaled_data = scaler.fit_transform(data)
+    return scaled_data, scaler
+
+
+def inverse_scale(scaled, scaler, n_features):
+    """Inverse transform using stored scaler."""
+    return scaler.inverse_transform(scaled.reshape(-1, n_features))
+
+
+def create_sequences(data, seq_length):
+    """Transform the time series into overlapping sequences."""
+    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)

visualize.py

@@ -0,0 +1,59 @@
+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()
+
+
+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.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_error_hist(errors):
+    plt.hist(errors, bins=50)
+    plt.title("Error Distribution")
+    plt.xlabel("Absolute Error")
+    plt.ylabel("Frequency")
+    plt.grid(True)
+    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_loss_curve(losses):
+    plt.plot(losses)
+    plt.title("Training Loss Curve")
+    plt.xlabel("Epoch")
+    plt.ylabel("MSE Loss")
+    plt.grid(True)
+    plt.tight_layout()
+    plt.show()