# ---------------------------------------------------------------------------
# traffic_stream.py
#
# Capstone project: warning a city's traffic control room fifteen minutes
# before a road jams, from a non-stop stream of sensor readings.
#
# CBSE Class XII, Artificial Intelligence (843), Part C - Capstone Project.
# SDG 11 (Sustainable Cities and Communities).
#
# Data Science Methodology, 10 steps, with BIG DATA ideas: the file is never
# loaded whole. It is read in chunks, summarised map-and-reduce style, cleaned
# chunk by chunk, and learned from in two big-data ways: an online model
# (partial_fit) that learns from every chunk in turn, and a Random Forest
# trained on a 10% sample. Both are compared with a simple rule.
#
# Libraries: pandas, numpy, matplotlib, scikit-learn, joblib.
# ---------------------------------------------------------------------------

# %% Step 1 - Business Understanding
print("=" * 64)
print("TRAFFIC JAM EARLY WARNING")
print("=" * 64)
print("Problem : by the time the control room sees a jam, it has already formed")
print("Goal    : warn 15 minutes ahead, so signals can be retimed and traffic")
print("          diverted while it still helps")
print()

# %% Step 2 - Analytic Approach
print("Analytic approach : predictive - will this road jam in 15 minutes? (yes/no)")
print("Jam means         : average speed below 20 km/h")
print()

# %% Step 3 - Data Requirements
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import SGDClassifier
from sklearn.ensemble import RandomForestClassifier
import joblib

FILE = "traffic_sensors.csv"
CHUNK = 20000                      # rows held in memory at a time
TEST_FROM = "2025-07-28"           # the last week is kept for testing
FEATURES = ["speed_kmph", "speed_drop_15min", "vehicles", "occupancy_pct", "rush_hour", "weekend"]

# %% Step 4 - Data Collection
print("--- Step 4: Data Collection ---")
print("File size :", round(os.path.getsize(FILE) / 1e6, 1), "MB")
rows = 0
chunks = 0
for chunk in pd.read_csv(FILE, chunksize=CHUNK):
    rows = rows + len(chunk)
    chunks = chunks + 1
print("Readings  :", rows)
print("Read in", chunks, "chunks of at most", CHUNK, "rows - never all at once")
print()

# %% Step 5 - Data Understanding
# MAP: summarise each chunk on its own. REDUCE: add the summaries together.
print("--- Step 5: Data Understanding ---")
speed_total = None
error_codes = 0
blank_counts = 0
jam_total = None
for chunk in pd.read_csv(FILE, chunksize=CHUNK):
    bad = (chunk["speed_kmph"] < 0) | (chunk["speed_kmph"] > 150)
    error_codes = error_codes + int(bad.sum())
    blank_counts = blank_counts + int(chunk["vehicles"].isnull().sum())
    good = chunk[~bad].copy()
    t = pd.to_datetime(good["time"])
    good["hour"] = t.dt.hour
    good["day"] = np.where(t.dt.dayofweek >= 5, "Weekend", "Weekday")
    part = good.groupby(["day", "hour"])["speed_kmph"].agg(["sum", "count"])           # map
    speed_total = part if speed_total is None else speed_total.add(part, fill_value=0)  # reduce
    jams = good.groupby(["sensor", "hour"])["jam_in_15_min"].agg(["sum", "count"])
    jam_total = jams if jam_total is None else jam_total.add(jams, fill_value=0)
avg_speed = (speed_total["sum"] / speed_total["count"]).unstack(0)
jam_share = (jam_total["sum"] / jam_total["count"]).unstack(1)
overall = jam_total["sum"].sum() / jam_total["count"].sum()
print("Speed error codes (-1 or 255) :", error_codes)
print("Blank vehicle counts          :", blank_counts)
print("Readings followed by a jam    :", round(100 * overall, 1), "%")
print("Slowest weekday hour          :", str(int(avg_speed["Weekday"].idxmin())) + ":00, average",
      round(avg_speed["Weekday"].min(), 1), "km/h")
print("Slowest weekend hour          :", str(int(avg_speed["Weekend"].idxmin())) + ":00, average",
      round(avg_speed["Weekend"].min(), 1), "km/h")
print("Sensor with the most jams     :", jam_share.mean(axis=1).idxmax())
print()

# Chart 1: average speed through the day.
plt.figure(figsize=(8, 4))
plt.plot(avg_speed.index, avg_speed["Weekday"], marker="o", label="Weekday", color="#e15759")
plt.plot(avg_speed.index, avg_speed["Weekend"], marker="o", label="Weekend", color="#4e79a7")
plt.axhline(20, color="gray", linestyle="--", label="Jam below 20 km/h")
plt.title("Average speed through the day, all sensors")
plt.xlabel("Hour of the day")
plt.ylabel("km/h")
plt.xticks(range(0, 24, 2))
plt.legend()
plt.grid(True, linestyle="--", alpha=0.4)
plt.tight_layout()
plt.savefig("chart1_hourly_speed.png")
plt.close()

# Chart 2: where and when jams follow.
plt.figure(figsize=(10, 5))
plt.imshow(jam_share.values, cmap="Reds", aspect="auto")
plt.colorbar(label="Share of readings followed by a jam")
plt.yticks(range(len(jam_share)), jam_share.index)
plt.xticks(range(0, 24, 2))
plt.xlabel("Hour of the day")
plt.ylabel("Sensor")
plt.title("Where and when jams happen")
plt.tight_layout()
plt.savefig("chart2_jam_heatmap.png")
plt.close()

# %% Step 6 - Data Preparation
# Each chunk is cleaned on its own. A chunk's first rows need the readings
# just before them (for the 15-minute change), so the last three readings of
# every sensor are CARRIED from one chunk to the next.
print("--- Step 6: Data Preparation ---")


def prepare(chunk, carry):
    chunk = chunk.copy()
    chunk["new"] = True
    both = pd.concat([carry, chunk])
    before = len(both)
    both = both.drop_duplicates(subset=["time", "sensor"])
    duplicates = before - len(both)
    bad = (both["speed_kmph"] < 0) | (both["speed_kmph"] > 150)
    both.loc[bad, "speed_kmph"] = np.nan
    both = both.sort_values(["sensor", "time"], kind="stable")
    # An error code or a blank count is replaced by that sensor's previous reading.
    both["speed_kmph"] = both.groupby("sensor")["speed_kmph"].ffill()
    both["vehicles"] = both.groupby("sensor")["vehicles"].ffill()
    both["speed_drop_15min"] = both.groupby("sensor")["speed_kmph"].shift(3) - both["speed_kmph"]
    t = pd.to_datetime(both["time"])
    both["weekend"] = (t.dt.dayofweek >= 5).astype(int)
    both["rush_hour"] = (t.dt.hour.isin([8, 9, 10, 17, 18, 19, 20]) & (both["weekend"] == 0)).astype(int)
    new_carry = both.groupby("sensor").tail(3).drop(columns="speed_drop_15min").assign(new=False)
    ready = both[both["new"] & both["speed_drop_15min"].notna()].sort_values(["time", "sensor"])
    return ready, new_carry, duplicates


def stream(file):
    """Yield one cleaned chunk at a time."""
    carry = pd.DataFrame()
    for chunk in pd.read_csv(file, chunksize=CHUNK):
        ready, carry, duplicates = prepare(chunk, carry)
        yield ready, duplicates


scaler = StandardScaler()
counts = np.zeros(2)
sample = []
duplicates_removed = 0
cleaned_rows = 0
for ready, dups in stream(FILE):
    duplicates_removed = duplicates_removed + dups
    cleaned_rows = cleaned_rows + len(ready)
    train = ready[ready["time"] < TEST_FROM]
    if len(train) > 0:
        scaler.partial_fit(train[FEATURES])                       # learns means bit by bit
        counts = counts + np.bincount(train["jam_in_15_min"], minlength=2)
        sample.append(train.sample(frac=0.1, random_state=0))    # a 10% sample for the forest
print("Readings sent twice, removed   :", duplicates_removed)
print("Error codes and blanks         : replaced by the sensor's previous reading")
print("Clean readings                 :", cleaned_rows)
print("Training readings (3 weeks)    :", int(counts.sum()))
print("  of which followed by a jam   :", int(counts[1]))
weight = counts[0] / counts[1]
print("Weight given to a jam reading  :", round(weight, 1), "(for the online model)")
print()

# Chart 3: one day at the sensor with the most error codes, before and after cleaning.
raw = pd.read_csv(FILE, nrows=CHUNK)
clean0, _, _ = prepare(raw, pd.DataFrame())
faulty = raw[(raw["speed_kmph"] < 0) | (raw["speed_kmph"] > 150)]
bad_sensor = faulty["sensor"].value_counts().index[0] if len(faulty) > 0 else raw["sensor"].iloc[0]
day = faulty[faulty["sensor"] == bad_sensor]["time"].iloc[0][:10] if len(faulty) > 0 else raw["time"].iloc[0][:10]
r = raw[(raw["sensor"] == bad_sensor) & raw["time"].str.startswith(day)].drop_duplicates(subset="time")
c = clean0[(clean0["sensor"] == bad_sensor) & clean0["time"].str.startswith(day)]
e = r[(r["speed_kmph"] < 0) | (r["speed_kmph"] > 150)]
plt.figure(figsize=(9, 4))
plt.plot(pd.to_datetime(r["time"]), r["speed_kmph"], color="#bab0ac", linewidth=3, label="As received")
plt.plot(pd.to_datetime(c["time"]), c["speed_kmph"], color="#4e79a7", label="Cleaned")
plt.scatter(pd.to_datetime(e["time"]), e["speed_kmph"], color="#e15759", zorder=3, label="Error code")
plt.title("Sensor " + bad_sensor + ", " + day + ": error codes replaced by the previous reading")
plt.ylabel("Speed, km/h")
plt.legend()
plt.tight_layout()
plt.savefig("chart3_cleaning.png")
plt.close()

# %% Step 7 - Modelling
print("--- Step 7: Modelling ---")
# average=True keeps a running average of the weights, so the online model
# does not swing with whichever chunk it saw last.
online = SGDClassifier(loss="log_loss", class_weight={0: 1, 1: weight}, average=True, random_state=0)
for epoch in range(3):                                            # three passes over the stream
    for ready, _ in stream(FILE):
        train = ready[ready["time"] < TEST_FROM]
        if len(train) > 0:
            online.partial_fit(scaler.transform(train[FEATURES]), train["jam_in_15_min"], classes=[0, 1])
print("Online model    : logistic regression learned with partial_fit, chunk by chunk, 3 passes")
small = pd.concat(sample)
forest = RandomForestClassifier(n_estimators=200, min_samples_leaf=5, class_weight="balanced", random_state=0)
forest.fit(small[FEATURES], small["jam_in_15_min"])
print("Random Forest   : trained on a 10% sample of", len(small), "readings held in memory")
print("Rule            : warn if the road is already jammed now")
print("Online model weights:")
for f, w in zip(FEATURES, online.coef_[0]):
    print("  ", f.ljust(18), round(w, 2))
print()

# %% Step 8 - Evaluation
print("--- Step 8: Evaluation ---")
names = ["Rule", "Online model", "Random Forest"]
tally = {n: np.zeros((2, 2), dtype=int) for n in names}          # [true][predicted]
onset = {n: np.zeros(3, dtype=int) for n in names}               # new jams warned, new jams, false alarms
by_hour = np.zeros((24, 2))                                       # jams caught, jams total
test_rows = 0
for ready, _ in stream(FILE):
    test = ready[ready["time"] >= TEST_FROM]
    if len(test) == 0:
        continue
    test_rows = test_rows + len(test)
    y = test["jam_in_15_min"].values
    clear = (test["speed_kmph"] >= 20).values                     # not jammed yet
    preds = {"Rule": (test["speed_kmph"] < 20).astype(int).values,
             "Online model": online.predict(scaler.transform(test[FEATURES])),
             "Random Forest": forest.predict(test[FEATURES])}
    for n in names:
        for a in (0, 1):
            for b in (0, 1):
                tally[n][a, b] += int(((y == a) & (preds[n] == b)).sum())
        onset[n] += [int((clear & (y == 1) & (preds[n] == 1)).sum()), int((clear & (y == 1)).sum()),
                     int((clear & (y == 0) & (preds[n] == 1)).sum())]
    hours = pd.to_datetime(test["time"]).dt.hour.values
    for h in range(24):
        jam = (hours == h) & (y == 1) & clear
        by_hour[h, 0] += (preds["Random Forest"][jam] == 1).sum()
        by_hour[h, 1] += jam.sum()
print("Test readings (last week) :", test_rows)
scores = {}
for n in names:
    tn, fp, fn, tp = tally[n].ravel()
    p = tp / (tp + fp)
    r = tp / (tp + fn)
    scores[n] = [p, r, 2 * p * r / (p + r), onset[n][0] / onset[n][1]]
    print(n.ljust(14), " TP", tp, " FN", fn, " FP", fp, " TN", tn)
print()
print("Model           precision  recall    F1")
for n in names:
    print(f"{n:<15} {scores[n][0]:9.3f} {scores[n][1]:7.3f} {scores[n][2]:6.3f}")
for n in names:
    print("Precision, recall, F1 -", n, ":", round(scores[n][0], 3), ",", round(scores[n][1], 3), ",", round(scores[n][2], 3))
print()
print("A warning matters most BEFORE a jam starts. Readings where the road")
print("was still clear but jammed 15 minutes later:", onset["Rule"][1])
for n in names:
    print(n + ", new jams warned of :", onset[n][0], "(" + str(round(100 * scores[n][3])) + "%)")
    print(n + ", false alarms on clear roads :", onset[n][2])
best = max(["Online model", "Random Forest"], key=lambda n: scores[n][2])
print("Best F1 :", best)
print()
rush = [8, 9, 10, 17, 18, 19, 20]
calm = [h for h in range(24) if h not in rush]
for label, hrs in [("in rush hours ", rush), ("at other times", calm)]:
    if by_hour[hrs, 1].sum() > 0:
        print("Random Forest, new jams warned of", label, ":", round(by_hour[hrs, 0].sum() / by_hour[hrs, 1].sum(), 3))
    else:
        print("Random Forest, new jams warned of", label, ": no new jams at these times")
print()

# Chart 4: the three approaches compared.
x = np.arange(3)
plt.figure(figsize=(9, 4))
labels = [("Precision", "#4e79a7"), ("Recall", "#e15759"), ("F1", "#59a14f"), ("New jams warned of", "#f28e2b")]
for k, (label, colour) in enumerate(labels):
    plt.bar(x + (k - 1.5) * 0.2, [scores[n][k] for n in names], width=0.2, label=label, color=colour)
plt.xticks(x, names)
plt.ylim(0, 1)
plt.title("Jam warnings on the test week")
plt.legend(loc="upper left", fontsize=8)
plt.tight_layout()
plt.savefig("chart4_models.png")
plt.close()

# Chart 5: which new jams the forest warns of, by hour.
plt.figure(figsize=(9, 4))
share = by_hour[:, 0] / np.maximum(by_hour[:, 1], 1)
plt.bar(range(24), share, color=["#e15759" if h in rush else "#9c9c9c" for h in range(24)])
plt.title("Share of new jams warned of in advance, by hour (red = weekday rush hours)")
plt.xlabel("Hour of the day")
plt.ylabel("Share warned of")
plt.xticks(range(0, 24, 2))
plt.ylim(0, 1)
plt.tight_layout()
plt.savefig("chart5_recall_by_hour.png")
plt.close()

# %% Step 9 - Deployment
# The saved model reads the stream as it arrives, and warns on the fly.
print("--- Step 9: Deployment ---")
joblib.dump(forest, "jam_model.pkl")
model = joblib.load("jam_model.pkl")
evening = None
for ready, _ in stream(FILE):
    part = ready[(ready["sensor"] == "J05") & (ready["time"] >= "2025-07-30 15:00") & (ready["time"] <= "2025-07-30 22:00")]
    if len(part) > 0:
        evening = part if evening is None else pd.concat([evening, part])
alerts = []
for i in range(len(evening)):
    reading = evening.iloc[[i]]                                   # one reading arrives
    if model.predict(reading[FEATURES])[0] == 1:
        alerts.append(reading["time"].iloc[0])
print("Sensor J05, 30 July, 3 pm to 10 pm :", len(evening), "readings streamed")
print("Warnings issued                    :", len(alerts))
print("First warning                      :", alerts[0] if alerts else "none")
jam_times = evening[evening["speed_kmph"] < 20]["time"]
print("Road first jammed                  :", jam_times.iloc[0] if len(jam_times) else "never")
print()

# Chart 6: the evening, as the control room would have seen it.
t = pd.to_datetime(evening["time"])
plt.figure(figsize=(9, 4))
plt.plot(t, evening["speed_kmph"], color="#4e79a7", label="Speed")
plt.axhline(20, color="gray", linestyle="--", label="Jam below 20 km/h")
warned = evening["time"].isin(alerts)
plt.scatter(t[warned], evening["speed_kmph"][warned], color="#f28e2b", zorder=3, label="Warning issued")
plt.title("Sensor J05, 30 July: warnings issued from the stream")
plt.ylabel("km/h")
plt.legend()
plt.tight_layout()
plt.savefig("chart6_stream_evening.png")
plt.close()

# %% Step 10 - Feedback
# Two festival shopping days, on a weekend: crowds when the model expects a
# quiet road. Day 1 is added to the training sample, and day 2 checked again.
print("--- Step 10: Feedback ---")
fest, _ = next(stream("festival_days.csv"))
day1 = fest[fest["time"] < "2025-10-19"]
day2 = fest[fest["time"] >= "2025-10-19"]


def check(m, d):
    p = m.predict(d[FEATURES])
    y = d["jam_in_15_min"].values
    clear = (d["speed_kmph"] >= 20).values
    return ((p == 1) & (y == 1)).sum() / (y == 1).sum(), ((p == 1) & (y == 1) & clear).sum() / ((y == 1) & clear).sum()


print("Festival day 1, new jams warned of :", round(check(model, day1)[1], 3))
r_before, n_before = check(model, day2)
print("Festival day 2 before retraining, recall             :", round(r_before, 3))
print("Festival day 2 before retraining, new jams warned of :", round(n_before, 3))
small2 = pd.concat([small, day1])
model2 = RandomForestClassifier(n_estimators=200, min_samples_leaf=5, class_weight="balanced", random_state=0)
model2.fit(small2[FEATURES], small2["jam_in_15_min"])
r_after, n_after = check(model2, day2)
print("Festival day 2 after adding day 1, recall             :", round(r_after, 3))
print("Festival day 2 after adding day 1, new jams warned of :", round(n_after, 3))

# Chart 7: festival days against an ordinary weekend, and the retraining.
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
fh = pd.to_datetime(fest["time"]).dt.hour
plt.plot(fest.groupby(fh)["speed_kmph"].mean(), color="#f28e2b", label="Festival weekend")
plt.plot(avg_speed["Weekend"], color="#4e79a7", label="Ordinary weekend")
plt.axhline(20, color="gray", linestyle="--")
plt.title("Average speed by hour")
plt.xlabel("Hour")
plt.ylabel("km/h")
plt.legend()
plt.subplot(1, 2, 2)
plt.bar(["Before retraining", "After adding day 1"], [n_before, n_after], color=["#9c9c9c", "#59a14f"])
plt.ylim(0, 1)
plt.title("Festival day 2: new jams warned of")
plt.tight_layout()
plt.savefig("chart7_festival.png")
plt.close()