Compare commits
2
Commits
No files matched your search
+1
-3
@@ -185,6 +185,4 @@ ipython_config.py
|
|||||||
|
|
||||||
|
|
||||||
.conda*/
|
.conda*/
|
||||||
|
dead/
|
||||||
data/
|
|
||||||
# for rn
|
|
||||||
File diff suppressed because it is too large.
Load diff
@@ -1,43 +0,0 @@
|
|||||||
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
|
|
||||||
+25
-15
@@ -1,12 +1,11 @@
|
|||||||
import urllib.parse
|
import urllib.parse
|
||||||
import urllib.request
|
import urllib.request
|
||||||
from datetime import datetime
|
from datetime import datetime
|
||||||
|
from global_config import DATA_FILE_NAME
|
||||||
data_file = "hanoi_weather_data.csv"
|
|
||||||
|
|
||||||
# Define parameters
|
# Define parameters
|
||||||
latitude, longitude = 21.04815477090267, 105.80123159165096 # lmao
|
latitude, longitude = 21.04815477090267, 105.80123159165096 # lmao
|
||||||
start_date = "1985-01-01"
|
start_date = "1981-01-01"
|
||||||
end_date = datetime.now().strftime("%Y-%m-%d") # Today's date
|
end_date = datetime.now().strftime("%Y-%m-%d") # Today's date
|
||||||
|
|
||||||
# Define desired parameters
|
# Define desired parameters
|
||||||
@@ -19,7 +18,19 @@ end_date = datetime.now().strftime("%Y-%m-%d") # Today's date
|
|||||||
# WD10M MERRA-2 Wind Direction at 10 Meters (Degrees)
|
# WD10M MERRA-2 Wind Direction at 10 Meters (Degrees)
|
||||||
# more at https://power.larc.nasa.gov/data-access-viewer/
|
# more at https://power.larc.nasa.gov/data-access-viewer/
|
||||||
weather_params = [
|
weather_params = [
|
||||||
|
# "T2M",
|
||||||
|
# "T2M_MAX",
|
||||||
|
# "T2M_MIN",
|
||||||
|
# "T2MDEW",
|
||||||
|
# "PRECTOTCORR",
|
||||||
"QV2M",
|
"QV2M",
|
||||||
|
# "RH2M",
|
||||||
|
# "PS",
|
||||||
|
# "CLOUD_AMT",
|
||||||
|
# "WS10M",
|
||||||
|
# "WD10M",
|
||||||
|
# "ALLSKY_KT",
|
||||||
|
# "ALLSKY_SFC_SW_DWN",
|
||||||
]
|
]
|
||||||
|
|
||||||
params = {
|
params = {
|
||||||
@@ -41,15 +52,14 @@ base_url = "https://power.larc.nasa.gov/api/temporal/daily/point"
|
|||||||
url = f"{base_url}?{query_string}"
|
url = f"{base_url}?{query_string}"
|
||||||
|
|
||||||
# Make the API request using urllib
|
# Make the API request using urllib
|
||||||
if __name__ == "__main__":
|
try:
|
||||||
try:
|
with urllib.request.urlopen(url) as response:
|
||||||
with urllib.request.urlopen(url) as response:
|
if response.status == 200:
|
||||||
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)
|
||||||
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:
|
||||||
else:
|
print(f"Failed to retrieve data: HTTP Status Code {response.status}")
|
||||||
print(f"Failed to retrieve data: HTTP Status Code {response.status}")
|
except urllib.error.URLError as e:
|
||||||
except urllib.error.URLError as e:
|
print(f"Failed to retrieve data: {e.reason}")
|
||||||
print(f"Failed to retrieve data: {e.reason}")
|
|
||||||
File diff suppressed because it is too large.
Load diff
@@ -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()
|
|
||||||
@@ -1 +0,0 @@
|
|||||||
DATA_FILE_NAME = "hanoi_weather_data.csv"
|
|
||||||
@@ -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()
|
|
||||||
@@ -0,0 +1 @@
|
|||||||
|
DATA_FILE_NAME = "data/hanoi_weather_data.csv"
|
||||||
File renamed without changes.
@@ -1,71 +0,0 @@
|
|||||||
# 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_original_vs_scaled,
|
|
||||||
plot_predictions,
|
|
||||||
plot_residuals,
|
|
||||||
plot_error_hist,
|
|
||||||
plot_scatter,
|
|
||||||
plot_loss_curve,
|
|
||||||
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)
|
|
||||||
|
|
||||||
# Step 2: Preprocess (scale and sequence)
|
|
||||||
scaled_data, scaler = scale_data(data)
|
|
||||||
plot_original_vs_scaled(data, scaled_data)
|
|
||||||
|
|
||||||
seq_length = 50
|
|
||||||
X, y = create_sequences(scaled_data, seq_length)
|
|
||||||
|
|
||||||
visualize_sequence_window(data["QV2M"], seq_length=10, start_idx=1000)
|
|
||||||
|
|
||||||
# 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
|
|
||||||
# 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,
|
|
||||||
n_epochs=10000,
|
|
||||||
)
|
|
||||||
|
|
||||||
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.to_numpy(), residuals)
|
|
||||||
plot_error_hist(errors)
|
|
||||||
plot_scatter(y_test.reshape(-1, 1), predictions.reshape(-1, 1))
|
|
||||||
plot_loss_curve(losses)
|
|
||||||
@@ -1,45 +0,0 @@
|
|||||||
import numpy as np
|
|
||||||
from sklearn.metrics import mean_squared_error, r2_score
|
|
||||||
|
|
||||||
|
|
||||||
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.
|
|
||||||
|
|
||||||
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):
|
|
||||||
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)
|
|
||||||
|
|
||||||
if epoch % 1000 == 0:
|
|
||||||
print(f"Epoch {epoch}: Loss = {loss:.6f}, LR: {lr:.6f}")
|
|
||||||
|
|
||||||
print(f"Final Loss: {losses[-1]:.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 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
|
|
||||||
@@ -1,26 +0,0 @@
|
|||||||
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)
|
|
||||||
@@ -4,40 +4,11 @@ import pandas as pd
|
|||||||
from sklearn.metrics import r2_score
|
from sklearn.metrics import r2_score
|
||||||
from sklearn.preprocessing import MinMaxScaler
|
from sklearn.preprocessing import MinMaxScaler
|
||||||
from sklearn.model_selection import train_test_split
|
from sklearn.model_selection import train_test_split
|
||||||
from sklearn.linear_model import LinearRegression
|
from sklearn.ensemble import RandomForestRegressor
|
||||||
from sklearn.metrics import mean_squared_error
|
from sklearn.metrics import mean_squared_error
|
||||||
import matplotlib.pyplot as plt
|
import matplotlib.pyplot as plt
|
||||||
# from global_config import DATA_FILE_NAME
|
from global_config import DATA_FILE_NAME
|
||||||
DATA_FILE_NAME = "hanoi_weather_data.csv"
|
from load_data import cut_headers
|
||||||
|
|
||||||
# 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
|
# 1. Load data
|
||||||
|
|
||||||
@@ -74,10 +45,16 @@ seq_length = 50
|
|||||||
X, y = create_sequences(scaled_data, seq_length)
|
X, y = create_sequences(scaled_data, seq_length)
|
||||||
|
|
||||||
# 4. Split into training and testing
|
# 4. Split into training and testing
|
||||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.02)
|
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.02, shuffle=False)
|
||||||
|
|
||||||
# 5. Train linear regression model
|
# 5. Train linear regression model
|
||||||
model = LinearRegression()
|
model = RandomForestRegressor(
|
||||||
|
n_estimators=200, # Number of trees
|
||||||
|
# max_features='sqrt', # Features considered at each split
|
||||||
|
# max_depth=30, # Maximum depth of each tree
|
||||||
|
# n_jobs=-1, # Use all CPU cores
|
||||||
|
random_state=42 # Reproducibility
|
||||||
|
)
|
||||||
model.fit(X_train, y_train)
|
model.fit(X_train, y_train)
|
||||||
|
|
||||||
# 6. Make predictions
|
# 6. Make predictions
|
||||||
@@ -96,15 +73,13 @@ predictions_rescaled = scaler.inverse_transform(predictions.reshape(-1, data.sha
|
|||||||
y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, data.shape[1]))
|
y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, data.shape[1]))
|
||||||
|
|
||||||
# 9. Plot results
|
# 9. Plot results
|
||||||
|
for i, column in enumerate(data.columns):
|
||||||
# Loss bar
|
plt.figure(figsize=(10, 5))
|
||||||
|
plt.plot(y_test_rescaled[:, i], label=f"Actual {column}")
|
||||||
|
plt.plot(predictions_rescaled[:, i], label=f"Predicted {column}")
|
||||||
# Prediction vs Actual
|
plt.title(f"Actual vs. Predicted for {column}")
|
||||||
plt.plot(y_test_rescaled, label="Actual")
|
plt.xlabel("Sample")
|
||||||
plt.plot(predictions_rescaled, label="Predicted")
|
plt.ylabel("Value")
|
||||||
plt.title("Actual vs. Predicted")
|
plt.legend()
|
||||||
plt.xlabel("Sample")
|
plt.tight_layout()
|
||||||
plt.ylabel("Value")
|
plt.show()
|
||||||
plt.legend()
|
|
||||||
plt.show()
|
|
||||||
-137
@@ -1,137 +0,0 @@
|
|||||||
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.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_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()
|
|
||||||
|
|
||||||
|
|
||||||
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 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).
|
|
||||||
|
|
||||||
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()
|
|
||||||
Reference in new issue
Block a user