1"""The three experiments (>= 5 seeds each, mean +- std, tables + plots).2 3 python run_experiments.py # all three, 10 seeds4 python run_experiments.py --seeds 5 --only 15 6Exp 1 A* vs UCS vs Greedy for 6..10 stops: nodes expanded/generated, cost, time7Exp 2 BT vs BT+FC vs BT+FC+AC-3 for 15/20/30 orders: assignments, backtracks, time8Exp 3 NN vs historical-average vs distance-based travel times: lateness in the replay9"""10import argparse11import json12import math13import shutil14import time15 16import numpy as np17import pandas as pd18 19from delivery import config as C20from delivery import plots as P21from delivery.csp import AssignmentCSP22from delivery.eval import ReplaySampler, replay23from delivery.instances import generate_instance, generate_route_instance24from delivery.integration import plan25from delivery.pipeline import load_all, time_model, zone_weights26from delivery.search import RouteProblem, brute_force, solve_route27 28RES = C.ROOT / "report" / "results"29FIG = C.ROOT / "report" / "figures"30 31 32def ms(x):33 x = np.asarray(x, float)34 return f"{x.mean():.4g} ± {x.std():.3g}"35 36 37def table(df, by, cols):38 rows = []39 for key, g in df.groupby(by, sort=False):40 key = key if isinstance(key, tuple) else (key,)41 rows.append(list(key) + [ms(g[c]) for c in cols])42 out = pd.DataFrame(rows, columns=list(by) + cols)43 return out44 45 46def md(df):47 cols = list(df.columns)48 lines = ["| " + " | ".join(cols) + " |", "|" + "---|" * len(cols)]49 lines += ["| " + " | ".join(str(v) for v in r) + " |" for r in df.values]50 return "\n".join(lines)51 52 53# ------------------------------------------------------------------------------ exp 154def exp1(mats, seeds):55 rows, bf_checks = [], []56 algos = [("UCS", "ucs", "zero"), ("Greedy", "greedy", "admissible"),57 ("A* (max-incoming h)", "astar", "max_in"), ("A* (sum/wait h)", "astar", "admissible")]58 for n in range(6, 11):59 for s in seeds:60 tm = time_model(mats, "nn", s % 5)61 inst = generate_route_instance(n, 1000 * n + s, tm, n_zones=tm.Z)62 prob = RouteProblem(inst.orders, inst.depot, inst.vans[0].shift_start, tm)63 res = {}64 for label, algo, h in algos:65 r = solve_route(prob, algo, heuristic=h)66 res[label] = r67 rows.append({"stops": n, "seed": s, "algorithm": label, "expanded": r.stats["expanded"],68 "generated": r.stats["generated"], "cost_min": r.cost,69 "time_ms": 1000 * r.stats["time_s"]})70 opt = res["UCS"].cost71 for row in rows[-len(algos):]:72 row["gap_pct"] = 100 * (row["cost_min"] - opt) / opt73 if n <= 8:74 t0 = time.perf_counter()75 bf = brute_force(prob)76 bf_checks.append({"stops": n, "seed": s, "brute_force": bf.cost,77 "astar": res["A* (sum/wait h)"].cost,78 "match": abs(bf.cost - res["A* (sum/wait h)"].cost) < 1e-6,79 "bf_time_ms": 1000 * (time.perf_counter() - t0),80 "permutations": bf.stats["permutations"]})81 df = pd.DataFrame(rows)82 df.to_csv(RES / "exp1_search.csv", index=False)83 bf = pd.DataFrame(bf_checks)84 bf.to_csv(RES / "exp1_bruteforce_check.csv", index=False)85 t = table(df, ["stops", "algorithm"], ["expanded", "generated", "cost_min", "gap_pct", "time_ms"])86 87 fig, axes = P.plt.subplots(1, 3, figsize=(14, 4))88 for k, (label, _, _) in enumerate(algos):89 g = df[df.algorithm == label].groupby("stops")90 for ax, col in zip(axes, ["expanded", "time_ms", "gap_pct"]):91 m, sd = g[col].mean(), g[col].std(ddof=0)92 ax.errorbar(m.index + (k - 1.5) * 0.06, m.values, yerr=sd.values, color=P.SERIES[k],93 marker="o", ms=5, capsize=3, label=label)94 axes[0].set_yscale("log"); axes[0].set_title("Nodes expanded (log)")95 axes[1].set_yscale("log"); axes[1].set_title("Run time (ms, log)")96 axes[2].set_title("Cost gap to optimum (%)")97 for ax in axes:98 ax.set_xlabel("stops per van"); ax.set_xticks(range(6, 11))99 axes[0].legend(loc="upper left")100 P.save(fig, C.OUTPUTS / "exp1_search.png")101 return t, bf102 103 104# ------------------------------------------------------------------------------ exp 2105def exp2(mats, zones, seeds):106 rows = []107 # "rush" instances: all windows 09:00-12:00, 15-45 min wide, capacity with zero slack.108 # Variables in input order for the three-way comparison (a fixed, uninformed order shows109 # the pruning itself); a 4th row adds the most-constrained-first order on top.110 solvers = [("BT", False, False, "input"), ("BT+FC", True, False, "input"),111 ("BT+FC+AC3", True, True, "input"), ("BT+FC+AC3 +degree order", True, True, "degree")]112 w = zone_weights(zones)113 for n in (15, 20, 30):114 for s in seeds:115 inst = generate_instance(n, s, n_zones=len(zones), zone_weights=w, cap_slack=1.0,116 max_stops=min(10, math.ceil(n / 3) + 1),117 day=(9 * 60, 12 * 60), widths=(15, 30, 45))118 tm = time_model(mats, "nn", inst.weekday)119 csp = AssignmentCSP(inst, tm)120 for label, fc, ac3, vo in solvers:121 r = csp.solve(fc=fc, ac3=ac3, max_assignments=200_000, time_limit=30, var_order=vo)122 rows.append({"orders": n, "seed": s, "solver": label, "status": r.status,123 "assignments": r.stats["assignments"], "backtracks": r.stats["backtracks"],124 "checks": r.stats["checks"], "fc_prunes": r.stats["fc_prunes"],125 "ac3_prunes": r.stats["ac3_prunes"], "time_ms": 1000 * r.stats["time_s"],126 "nn_edges": len(csp.edges)})127 df = pd.DataFrame(rows)128 df.to_csv(RES / "exp2_csp.csv", index=False)129 t = table(df, ["orders", "solver"], ["assignments", "backtracks", "checks", "time_ms", "nn_edges"])130 t.insert(3, "assignments median", df.groupby(["orders", "solver"], sort=False).assignments.median().values)131 st = df.groupby(["orders", "solver"], sort=False).status132 t["solved/unsat/limit"] = st.apply(lambda s: f"{(s == 'solved').sum()}/{(s == 'unsat').sum()}/{(s == 'limit').sum()}").values133 134 fig, axes = P.plt.subplots(1, 3, figsize=(14, 4))135 x = np.arange(3)136 for k, (label, _, _, _) in enumerate(solvers):137 g = df[df.solver == label].groupby("orders", sort=True)138 for ax, col in zip(axes, ["assignments", "backtracks", "time_ms"]):139 m = g[col].median()140 q1, q3 = g[col].quantile(0.25), g[col].quantile(0.75)141 ax.bar(x + (k - 1.5) * 0.2, m.values + (1 if col != "time_ms" else 0), 0.19,142 yerr=[m.values - q1.values, q3.values - m.values],143 color=P.SERIES[k], label=label, capsize=3, edgecolor=P.SURFACE, linewidth=2,144 error_kw={"ecolor": P.MUTED, "elinewidth": 1})145 titles = ["Assignments (median, IQR, log)", "Backtracks (+1, median, log)", "Time (ms, median, log)"]146 for ax, ti in zip(axes, titles):147 ax.set_yscale("log"); ax.set_title(ti); ax.set_xticks(x); ax.set_xticklabels(["15", "20", "30"])148 ax.set_xlabel("orders")149 axes[0].legend(loc="upper left")150 P.save(fig, C.OUTPUTS / "exp2_csp.png")151 return t152 153 154# ------------------------------------------------------------------------------ exp 3155def exp3(mats, zones, trips, seeds, reps=30):156 sampler = ReplaySampler(trips[trips.split == "test"])157 w = zone_weights(zones)158 rows = []159 scenarios = {"standard windows 30-120 min": (30, 45, 60, 90, 120), "tight windows 15-30 min": (15, 30)}160 for scen, widths in scenarios.items():161 for s in seeds:162 inst = generate_instance(20, 100 + s, n_zones=len(zones), zone_weights=w, widths=widths)163 for kind in ("nn", "hist", "dist"):164 tm = time_model(mats, kind, inst.weekday)165 res = plan(inst, tm, suggest=False)166 rp = replay(inst, res, sampler, tm, n_reps=reps, seed=s)167 planned_time = sum(r.finish - inst.vans[v].shift_start for v, r in res.routes.items() if r.sequence)168 rows.append({"scenario": scen, "seed": s, "model": tm.name, "plan_status": res.status, "loops": len(res.loops),169 "learned_rules": len(res.learned_edges), "planned_late": res.n_late,170 "replay_late": rp["late_mean"], "replay_late_pct": 100 * rp["late_frac"],171 "replay_late_minutes": rp["late_minutes_mean"], "unserved": rp["unserved"],172 "planned_route_min": planned_time})173 df = pd.DataFrame(rows)174 df.to_csv(RES / "exp3_lateness.csv", index=False)175 t = table(df, ["scenario", "model"], ["replay_late", "replay_late_pct", "replay_late_minutes", "planned_late",176 "planned_route_min", "learned_rules"])177 t["plans ok"] = df.groupby(["scenario", "model"], sort=False).plan_status.apply(178 lambda s: f"{(s == 'ok').sum()}/{len(s)}").values179 180 fig, axes = P.plt.subplots(1, 2, figsize=(12, 4.2), sharey=False)181 models = list(df.model.unique())182 for a, (ax, scen) in enumerate(zip(axes, scenarios)):183 d = df[df.scenario == scen]184 for k, m in enumerate(models):185 g = d[d.model == m]186 ax.bar(k, g.replay_late_pct.mean(), 0.6, color=P.SERIES[k], edgecolor=P.SURFACE, linewidth=2)187 ax.scatter(np.full(len(g), k) + np.linspace(-0.18, 0.18, len(g)), g.replay_late_pct, s=12,188 color=P.INK, zorder=4)189 ax.annotate(f"{g.replay_late_pct.mean():.1f}%", (k, g.replay_late_pct.mean()), xytext=(14, 4),190 textcoords="offset points", color=P.INK2)191 ax.set_xticks(range(len(models))); ax.set_xticklabels(models)192 ax.set_title(f"Late in replay: {scen}")193 ax.set_ylabel("% of orders late (bar = mean, dots = seeds)")194 P.save(fig, C.OUTPUTS / "exp3_lateness.png")195 return t, sampler.backoffs196 197 198def main():199 ap = argparse.ArgumentParser()200 ap.add_argument("--seeds", type=int, default=10)201 ap.add_argument("--only", type=int, default=0)202 args = ap.parse_args()203 seeds = list(range(args.seeds))204 RES.mkdir(parents=True, exist_ok=True); FIG.mkdir(parents=True, exist_ok=True); C.OUTPUTS.mkdir(exist_ok=True)205 trips, zones, dist, mats = load_all()206 summary = {}207 if args.only in (0, 1):208 print("Experiment 1: search"); t0 = time.time()209 t, bf = exp1(mats, seeds)210 print(t.to_string(index=False)); print(f"A* == brute force: {int(bf.match.sum())}/{len(bf)} ({time.time()-t0:.0f}s)")211 summary["exp1"] = md(t) + f"\n\nA* vs brute force (6-8 stops): {int(bf.match.sum())}/{len(bf)} identical optimal costs.\n"212 if args.only in (0, 2):213 print("Experiment 2: CSP"); t0 = time.time()214 t = exp2(mats, zones, seeds)215 print(t.to_string(index=False)); print(f"({time.time()-t0:.0f}s)")216 summary["exp2"] = md(t)217 if args.only in (0, 3):218 print("Experiment 3: lateness replay"); t0 = time.time()219 t, back = exp3(mats, zones, trips, seeds)220 print(t.to_string(index=False)); print("replay back-offs:", dict(back), f"({time.time()-t0:.0f}s)")221 summary["exp3"] = md(t) + f"\n\nReplay pool back-offs: {dict(back)}\n"222 for name, txt in summary.items():223 (RES / f"{name}_table.md").write_text(txt, encoding="utf-8")224 for f in C.OUTPUTS.glob("exp*.png"):225 shutil.copy(f, FIG / f.name)226 print("tables -> report/results/, figures -> outputs/ and report/figures/")227 228 229if __name__ == "__main__":230 main()