# ---------------------------------------------------------------------------
# power_demand.py
#
# Capstone project: forecasting tomorrow's peak electricity demand for a
# city, so the distribution company can buy and schedule the right amount
# of power - no blackouts, and no expensive power bought and wasted.
#
# CBSE Class XII, Artificial Intelligence (843), Part C - Capstone Project.
# SDG 7 (Affordable and Clean Energy) and SDG 13 (Climate Action).
#
# Data Science Methodology, 10 steps. REGRESSION on a time series. Three
# models are compared - Linear Regression, a Random Forest and a neural
# network (a multi-layer perceptron) - and, because this is a forecast,
# they are trained on the past and tested on a later year, then checked
# with time-series cross validation.
#
# Libraries: pandas, numpy, matplotlib, scikit-learn, joblib.
# ---------------------------------------------------------------------------

# %% Step 1 - Business Understanding
print("=" * 64)
print("TOMORROW'S ELECTRICITY DEMAND")
print("=" * 64)
print("Problem : power for tomorrow is bought today; a bad guess means")
print("          load-shedding or paying for power nobody uses")
print("Goal    : forecast tomorrow's peak demand (MW) from weather and calendar")
print()

# %% Step 2 - Analytic Approach
print("Analytic approach : predictive - regression on a time series")
print()

# %% Step 3 - Data Requirements
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import joblib

pd.set_option("display.width", 130)
pd.set_option("display.max_columns", 20)

# %% Step 4 - Data Collection
data = pd.read_csv("power_demand.csv", parse_dates=["date"])
print("Rows collected :", len(data))
print("Dates          :", data["date"].min().date(), "to", data["date"].max().date())
print(data.head(6))
print()

# %% Step 5 - Data Understanding
print("--- Step 5: Data Understanding ---")
print("Missing values:")
print(data.isnull().sum().to_string())
print("Days with demand recorded as 0 :", (data["demand_mw"] == 0).sum())
print()
print(data.describe().round(1))
print()

ok = data[data["demand_mw"] > 0]

# Chart 1: demand every day for three years.
plt.figure(figsize=(10, 4))
plt.plot(ok["date"], ok["demand_mw"], color="#f28e2b", linewidth=0.8)
plt.title("Daily peak demand, 2022-2024")
plt.xlabel("Date")
plt.ylabel("Demand (MW)")
plt.grid(True, linestyle="--", alpha=0.3)
plt.tight_layout()
plt.savefig("chart1_demand.png")
plt.close()

# Chart 2: temperature against demand.
plt.figure(figsize=(7, 5))
plt.scatter(ok["max_temp_c"], ok["demand_mw"], s=8, alpha=0.5, color="#e15759")
plt.title("Maximum temperature and demand")
plt.xlabel("Maximum temperature (deg C)")
plt.ylabel("Demand (MW)")
plt.grid(True, linestyle="--", alpha=0.3)
plt.tight_layout()
plt.savefig("chart2_temperature.png")
plt.close()

# Chart 3: average demand on each day of the week.
weekday = ok.groupby(ok["date"].dt.dayofweek)["demand_mw"].mean()
names = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"]
print("Average demand by day of the week (MW):")
for i in range(7):
    print("  " + names[i], round(weekday[i]))
print()
plt.figure(figsize=(7, 4))
plt.bar(names, weekday.values, color=["#4e79a7"] * 5 + ["#9c9c9c", "#9c9c9c"])
plt.ylim(weekday.min() - 300, weekday.max() + 100)
plt.title("Average demand by day of the week")
plt.ylabel("Demand (MW)")
plt.tight_layout()
plt.savefig("chart3_weekday.png")
plt.close()

# %% Step 6 - Data Preparation
print("--- Step 6: Data Preparation ---")
before = len(data)
data = data.drop_duplicates(subset="date").sort_values("date").reset_index(drop=True)
print("Duplicate dates removed        :", before - len(data))

# A demand of 0 is a grid failure, not a real demand: it becomes missing.
# Missing days are then filled by interpolating between the days either
# side, which suits a smooth day-to-day series.
data["demand_mw"] = pd.to_numeric(data["demand_mw"]).replace(0, np.nan)
gaps = data["demand_mw"].isnull().sum()
data["demand_mw"] = data["demand_mw"].interpolate()
print("Missing or zero days filled    :", gaps, "(interpolated)")

# Feature engineering - all from information known the day before.
data["day_of_week"] = data["date"].dt.dayofweek
data["is_sunday"] = (data["day_of_week"] == 6).astype(int)
data["is_saturday"] = (data["day_of_week"] == 5).astype(int)
data["month"] = data["date"].dt.month
data["cooling_degrees"] = np.maximum(0, data["max_temp_c"] - 27)       # heat above 27 drives ACs
data["demand_yesterday"] = data["demand_mw"].shift(1)
data["demand_last_week"] = data["demand_mw"].shift(7)
data = data.dropna().reset_index(drop=True)
features = ["max_temp_c", "cooling_degrees", "humidity_pct", "holiday", "is_saturday", "is_sunday",
            "month", "demand_yesterday", "demand_last_week"]
print("Rows after making lag features :", len(data), "(first 7 days have no 'last week')")
print("Features used                  :", features)
print()

# Chart 4: demand in each month, as box plots.
plt.figure(figsize=(9, 4.3))
plt.boxplot([data[data["month"] == m]["demand_mw"] for m in range(1, 13)],
            labels=["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"])
plt.title("Demand in each month (2022-2024)")
plt.ylabel("Demand (MW)")
plt.grid(True, axis="y", linestyle="--", alpha=0.3)
plt.tight_layout()
plt.savefig("chart4_months.png")
plt.close()

# %% Step 7 - Modelling
print("--- Step 7: Modelling ---")
# A forecast must be tested on a future the model has not seen, so the
# rows are NOT shuffled: train on 2022-2023, test on 2024.
train = data[data["date"].dt.year <= 2023]
test = data[data["date"].dt.year == 2024]
X_train, y_train = train[features], train["demand_mw"]
X_test, y_test = test[features], test["demand_mw"]
print("Training days (2022-2023) :", len(train))
print("Testing days (2024)       :", len(test))
print()

models = {
    "Linear Regression": LinearRegression(),
    "Random Forest": RandomForestRegressor(n_estimators=300, min_samples_leaf=2, random_state=0),
    "Neural Network": make_pipeline(StandardScaler(),
                                    MLPRegressor(hidden_layer_sizes=(32, 16), max_iter=3000,
                                                 learning_rate_init=0.01, random_state=0)),
}
for name in models:
    models[name].fit(X_train, y_train)
print("Neural network layers : input 9 -> 32 -> 16 -> output 1")
print()

# %% Step 8 - Evaluation
print("--- Step 8: Evaluation ---")
average = y_test.mean()
predictions = {}
for name in models:
    pred = models[name].predict(X_test)
    predictions[name] = pred
    rmse = np.sqrt(mean_squared_error(y_test, pred))
    print(name)
    print("  " + name + " MAE  :", round(mean_absolute_error(y_test, pred), 1), "MW")
    print("  " + name + " RMSE :", round(rmse, 1), "MW")
    print("  " + name + " R2   :", round(r2_score(y_test, pred), 3))
    print("  " + name + " RMSE as % of average demand :", round(rmse / average * 100, 2), "%")
naive = test["demand_yesterday"]
print("Baseline 'same as yesterday' RMSE :", round(np.sqrt(mean_squared_error(y_test, naive)), 1), "MW")
print()

# Time-series cross validation: each fold trains on everything before it
# and tests on the block that follows - never on the past.
tscv = TimeSeriesSplit(n_splits=5)
cv_rmse = {name: [] for name in models}
X_all, y_all = data[features], data["demand_mw"]
for train_idx, test_idx in tscv.split(X_all):
    for name in models:
        m = {"Linear Regression": LinearRegression(),
             "Random Forest": RandomForestRegressor(n_estimators=300, min_samples_leaf=2, random_state=0),
             "Neural Network": make_pipeline(StandardScaler(),
                                             MLPRegressor(hidden_layer_sizes=(32, 16), max_iter=3000,
                                                          learning_rate_init=0.01, random_state=0))}[name]
        m.fit(X_all.iloc[train_idx], y_all.iloc[train_idx])
        p = m.predict(X_all.iloc[test_idx])
        cv_rmse[name].append(np.sqrt(mean_squared_error(y_all.iloc[test_idx], p)))
print("Time-series cross validation, RMSE (MW) in each fold:")
for name in models:
    print("  " + name.ljust(18), np.round(cv_rmse[name], 1), " mean", round(np.mean(cv_rmse[name]), 1))
print()

# Chart 5: actual and predicted demand through 2024.
plt.figure(figsize=(10, 4.3))
plt.plot(test["date"], y_test, color="#999999", linewidth=1, label="Actual")
plt.plot(test["date"], predictions["Linear Regression"], color="#4e79a7", linewidth=1, label="Linear Regression")
plt.title("2024: actual demand and the Linear Regression forecast")
plt.xlabel("Date")
plt.ylabel("Demand (MW)")
plt.legend()
plt.grid(True, linestyle="--", alpha=0.3)
plt.tight_layout()
plt.savefig("chart5_forecast_2024.png")
plt.close()

# Chart 6: error in each time-series fold for the three models.
x = np.arange(5)
plt.figure(figsize=(8, 4.3))
plt.bar(x - 0.27, cv_rmse["Linear Regression"], width=0.27, label="Linear Regression", color="#4e79a7")
plt.bar(x, cv_rmse["Random Forest"], width=0.27, label="Random Forest", color="#59a14f")
plt.bar(x + 0.27, cv_rmse["Neural Network"], width=0.27, label="Neural Network", color="#b07aa1")
plt.xticks(x, ["Fold 1", "Fold 2", "Fold 3", "Fold 4", "Fold 5"])
plt.ylabel("RMSE (MW)")
plt.title("Time-series cross validation")
plt.legend()
plt.tight_layout()
plt.savefig("chart6_tscv.png")
plt.close()

# Chart 7: what the linear model's coefficients mean, per unit.
coef = pd.Series(models["Linear Regression"].coef_, index=features)
print("Linear Regression coefficients:")
print(coef.round(2).to_string())
print()
summer = test[(test["date"].dt.month >= 4) & (test["date"].dt.month <= 6)]
plt.figure(figsize=(10, 4.3))
plt.plot(summer["date"], summer["demand_mw"], color="#999999", label="Actual")
for name, col in [["Linear Regression", "#4e79a7"], ["Random Forest", "#59a14f"], ["Neural Network", "#b07aa1"]]:
    plt.plot(summer["date"], models[name].predict(summer[features]), color=col, linewidth=1, label=name)
plt.title("Summer 2024 close up: April to June")
plt.ylabel("Demand (MW)")
plt.legend()
plt.grid(True, linestyle="--", alpha=0.3)
plt.tight_layout()
plt.savefig("chart7_summer.png")
plt.close()

# %% Step 9 - Deployment
print("--- Step 9: Deployment ---")
best = models["Linear Regression"]
joblib.dump(best, "demand_model.pkl")
model = joblib.load("demand_model.pkl")
last = data.iloc[-1]
last_week = data.iloc[-7]
tomorrow = pd.DataFrame({"max_temp_c": [24.0], "cooling_degrees": [0.0], "humidity_pct": [58], "holiday": [0],
                         "is_saturday": [0], "is_sunday": [0], "month": [1],
                         "demand_yesterday": [last["demand_mw"]], "demand_last_week": [last_week["demand_mw"]]})
print("Forecast for 1 January 2025 (24 deg C forecast) :", round(model.predict(tomorrow)[0]), "MW")
print()

# %% Step 10 - Feedback
print("--- Step 10: Feedback ---")
# Each day of January 2025 is forecast the evening before, then compared
# with what the meters actually recorded.
jan = pd.read_csv("january_2025.csv", parse_dates=["date"])
history = pd.concat([data[["date", "demand_mw"]], jan[["date", "demand_mw"]]]).reset_index(drop=True)
jan["day_of_week"] = jan["date"].dt.dayofweek
jan["is_sunday"] = (jan["day_of_week"] == 6).astype(int)
jan["is_saturday"] = (jan["day_of_week"] == 5).astype(int)
jan["month"] = 1
jan["cooling_degrees"] = np.maximum(0, jan["max_temp_c"] - 27)
jan["demand_yesterday"] = history["demand_mw"].shift(1).iloc[-len(jan):].values
jan["demand_last_week"] = history["demand_mw"].shift(7).iloc[-len(jan):].values
jan_pred = model.predict(jan[features])
jan_rmse = np.sqrt(mean_squared_error(jan["demand_mw"], jan_pred))
test_rmse = np.sqrt(mean_squared_error(y_test, predictions["Linear Regression"]))
jan_2024 = test[test["date"].dt.month == 1]
jan_2024_rmse = np.sqrt(mean_squared_error(jan_2024["demand_mw"], model.predict(jan_2024[features])))
print("Days forecast in January 2025     :", len(jan))
print("RMSE, January 2025                :", round(jan_rmse, 1), "MW")
print("RMSE, January 2024 (test year)    :", round(jan_2024_rmse, 1), "MW")
print("RMSE, all of 2024                 :", round(test_rmse, 1), "MW")
print("Average error in January 2025     :", round((jan["demand_mw"] - jan_pred).mean(), 1),
      "MW (actual minus forecast)")
if jan_rmse > 1.25 * jan_2024_rmse:
    print("Feedback: January 2025 errors are well above last January's - retrain.")
else:
    print("Feedback: January 2025 is in line with last January - keep the model, retrain each quarter.")

# Chart 8: January 2025, forecast against actual.
plt.figure(figsize=(9, 4))
plt.plot(jan["date"], jan["demand_mw"], marker="o", markersize=3, color="#999999", label="Actual")
plt.plot(jan["date"], jan_pred, marker="s", markersize=3, color="#4e79a7", label="Forecast the day before")
plt.title("Feedback: January 2025, day-ahead forecasts against the meters")
plt.ylabel("Demand (MW)")
plt.legend()
plt.grid(True, linestyle="--", alpha=0.3)
plt.tight_layout()
plt.savefig("chart8_january_2025.png")
plt.close()