# ---------------------------------------------------------------------------
# home_energy.py
#
# AI project on Statistical Data: predicting how many units of electricity a
# home will use in a summer month, from what is in the home.
#
# CBSE Class X, Artificial Intelligence (417), Part C project work.
# SDG 7 (Affordable and Clean Energy) and SDG 13 (Climate Action).
#
# This is a REGRESSION project: the answer is a number, not a group. Three
# models are compared, each a little cleverer than the last:
#   A. always guess the average home
#   B. a line of best fit through one feature (AC hours a day)
#   C. a weighted sum of every feature - each feature times a weight, plus a
#      bias, the same sum a perceptron works out in Unit 2
#
# They are judged the way the curriculum judges a regression model:
#   error = |actual - predicted|,  error rate = error / actual,
#   accuracy = 1 - error rate,  model accuracy = mean of all the accuracies.
#
# Libraries: pandas (the data), numpy (the arithmetic), matplotlib (charts).
# ---------------------------------------------------------------------------

# %% Stage 1 - Problem scoping
print("=" * 60)
print("HOW MANY UNITS WILL WE USE?")
print("=" * 60)
print("Who   : families paying the electricity bill")
print("What  : the summer bill arrives as a shock, too late to act")
print("Where : homes in a hot city, in the months the AC runs")
print("Why   : knowing in advance helps a family save power")
print()

# %% Stage 2 - Data acquisition
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

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

data = pd.read_csv("home_energy.csv")

print("--- Stage 2: Data acquisition ---")
print("Rows read    :", len(data))
print("Columns read :", len(data.columns))
print()
print("The first 10 rows:")
print(data.head(10))
print()

# %% Stage 3 - Data exploration
print("--- Stage 3: Data exploration ---")

# A home with no meter reading cannot teach the model anything.
missing = data["units"].isnull().sum()
data = data.dropna(subset=["units"])
print("Rows with no meter reading       :", missing, "(removed)")

# A day has 24 hours, so 60 or 80 AC hours a day is a typing mistake.
impossible = (data["ac_hours"] > 24).sum()
data = data[data["ac_hours"] <= 24]
print("Rows with more than 24 AC hours  :", impossible, "(removed)")
print("Rows left                        :", len(data))
print()

# One new feature, made from two: how many AC-hours the home runs in a day.
# Two ACs for 5 hours use about as much power as one AC for 10.
data["ac_total_hours"] = data["ac_count"] * data["ac_hours"]

units = data["units"]
print("Units in a month - mean   :", round(np.mean(units), 1))
print("Units in a month - median :", np.median(units))
print("Units in a month - lowest :", units.min(), " highest :", units.max())
print()

by_members = data.groupby("members")["units"].mean().round(1)
print("Average units, by number of people in the home:")
print(by_members)
print()

# Chart 1: AC hours a day against units.
plt.figure(figsize=(7, 5))
plt.scatter(data["ac_total_hours"], data["units"], alpha=0.6, s=18, color="#f28e2b")
plt.title("AC hours a day and units used in the month")
plt.xlabel("AC hours a day (all ACs added up)")
plt.ylabel("Units (kWh) in the month")
plt.grid(True, linestyle="--", alpha=0.4)
plt.tight_layout()
plt.savefig("chart1_ac_hours.png")
plt.close()

# Chart 2: average units by family size.
plt.figure(figsize=(7, 4))
plt.bar(by_members.index, by_members.values, color="#4e79a7")
plt.title("Average units used, by number of people")
plt.xlabel("People in the home")
plt.ylabel("Average units (kWh)")
plt.tight_layout()
plt.savefig("chart2_members.png")
plt.close()

# %% Stage 4 - Modelling
print("--- Stage 4: Modelling ---")

features = ["members", "rooms", "ac_total_hours", "geyser"]
label = "units"

data = data.sample(frac=1, random_state=5)
cut = int(len(data) * 0.8)
train = data[:cut]
test = data[cut:]
print("Training rows :", len(train))
print("Testing rows  :", len(test))
print()

# Model A: always guess the average of the training homes.
average = train[label].mean()
predict_a = np.full(len(test), average)
print("Model A guesses every home will use", round(average, 1), "units")

# Model B: the straight line that fits AC hours and units best.
# np.polyfit with 1 finds the slope and the starting point of that line.
slope, start = np.polyfit(train["ac_total_hours"], train[label], 1)
predict_b = slope * test["ac_total_hours"] + start
print("Model B line : units =", round(start, 1), "+", round(slope, 1), "x AC hours a day")

# Model C: a weighted sum of all four features, plus a bias.
#     units = bias + w1 x members + w2 x rooms + w3 x AC hours + w4 x geyser
# np.linalg.lstsq finds the weights that make the predictions on the
# training homes as close as possible to the real units. The column of 1s
# is what lets it find the bias, exactly like the bias input of a perceptron.
X_train = train[features].copy()
X_train["bias"] = 1
weights = np.linalg.lstsq(X_train.values, train[label].values.astype(float), rcond=None)[0]

print()
print("What model C learned (one weight for each feature):")
for i in range(len(features)):
    print("  " + features[i].ljust(16), round(weights[i], 1))
print("  " + "bias".ljust(16), round(weights[4], 1))
print()

X_test = test[features].copy()
X_test["bias"] = 1
predict_c = X_test.values.dot(weights)

# %% Stage 5 - Evaluation
print("--- Stage 5: Evaluation ---")
actual = test[label].values.astype(float)

names = ["A: average", "B: one line", "C: weighted sum"]
model_accuracy = []
mean_error = []
for predicted in [predict_a, predict_b, predict_c]:
    error = np.abs(actual - predicted)          # error = |actual - predicted|
    error_rate = error / actual                 # error rate = error / actual
    accuracy = 1 - error_rate                   # accuracy = 1 - error rate
    model_accuracy.append(np.mean(accuracy) * 100)
    mean_error.append(np.mean(error))

for i in range(3):
    print("Model " + names[i].ljust(16), "average error :", str(round(mean_error[i], 1)).rjust(6),
          "units   accuracy :", round(model_accuracy[i], 1), "%")
print()
print("Average error of model C :", round(mean_error[2], 1), "units")
print("Accuracy of model C      :", round(model_accuracy[2], 1), "%")
print("Accuracy of model B      :", round(model_accuracy[1], 1), "%")
print("Accuracy of model A      :", round(model_accuracy[0], 1), "%")
print()

# Chart 3: model C's predictions against the real units.
plt.figure(figsize=(6, 6))
plt.scatter(actual, predict_c, alpha=0.7, s=20, color="#59a14f")
top = max(actual.max(), predict_c.max())
plt.plot([0, top], [0, top], color="#999999", linestyle="--", label="Perfect prediction")
plt.title("Model C: predicted against actual units")
plt.xlabel("Actual units")
plt.ylabel("Predicted units")
plt.legend()
plt.grid(True, linestyle="--", alpha=0.4)
plt.tight_layout()
plt.savefig("chart3_predicted_actual.png")
plt.close()

# Chart 4: the three models compared.
plt.figure(figsize=(7, 4))
plt.bar(names, model_accuracy, color=["#bab0ac", "#f28e2b", "#59a14f"])
plt.ylim(0, 100)
plt.title("Accuracy of the three models (testing homes)")
plt.ylabel("Model accuracy (%)")
for i in range(3):
    plt.text(i, model_accuracy[i] + 2, str(round(model_accuracy[i], 1)) + "%", ha="center")
plt.tight_layout()
plt.savefig("chart4_models.png")
plt.close()

# Chart 5: how far off each testing home was, for model C.
plt.figure(figsize=(7, 4))
plt.hist(actual - predict_c, bins=15, color="#59a14f")
plt.title("Model C: how far off each prediction was")
plt.xlabel("Actual minus predicted (units)")
plt.ylabel("Number of homes")
plt.tight_layout()
plt.savefig("chart5_errors.png")
plt.close()

# %% Stage 6 - Deployment
print("--- Stage 6: Deployment ---")
homes = pd.DataFrame({
    "members":        [4, 2, 6],
    "rooms":          [3, 1, 4],
    "ac_total_hours": [16, 0, 24],
    "geyser":         [1, 0, 1],
})
homes["bias"] = 1
forecast = homes.values.dot(weights)
for i in range(len(homes)):
    h = homes.iloc[i]
    print("People", h["members"], " rooms", h["rooms"], " AC hours a day", h["ac_total_hours"],
          " geyser", h["geyser"], "->", round(forecast[i]), "units")

# What would two hours less AC a day save the first home?
fewer = homes.iloc[[0]].copy()
fewer["ac_total_hours"] = fewer["ac_total_hours"] - 2
saving = forecast[0] - fewer.values.dot(weights)[0]
print("Saving from two AC hours a day less :", round(saving), "units a month for the first home")