1"""End-to-end demo: NN travel times -> CSP van assignment -> A* stop order, with the feedback loop.2 3 python run_demo.py # default dispatcher day (seed 4), figures in outputs/4 python run_demo.py --impossible # same day with one unreachable order -> "no solution" report5 python run_demo.py --seed 36 # another generated day6 7Prints the input, the loop log, the routes, a search trace, CSP pruning counters and8NN-on vs NN-off replay lateness; saves demo_*.png and demo_report.html to outputs/.9"""10import argparse11import html12import io13import math14from contextlib import redirect_stdout15 16import numpy as np17 18from delivery import config as C19from delivery import plots as P20from delivery.csp import AssignmentCSP, fmt_t21from delivery.data import read_zone_shapes22from delivery.eval import ReplaySampler, replay23from delivery.instances import generate_instance24from delivery.integration import describe_plan, plan, van_problem25from delivery.pipeline import load_all, time_model, zone_weights26from delivery.search import solve_route27 28DEPOT_ZONE = 7 # Midtown Center29plt = P.plt30 31 32# ------------------------------------------------------------------------------ figures33def stop_positions(inst, zones):34 """Centroid per stop, with stops in the same zone fanned out on a small circle."""35 xy = zones.sort_values("idx")[["cx_km", "cy_km"]].to_numpy()36 pos, seen = {}, {}37 for o in inst.orders:38 k = seen.get(o.zone, 0); seen[o.zone] = k + 139 ang = 2.4 * k40 r = 0.0 if k == 0 else 0.2841 pos[o.id] = xy[o.zone] + r * np.array([np.cos(ang), np.sin(ang)])42 return xy, pos43 44 45def draw_map(inst, res, zones, path):46 shapes_path = C.DATA_RAW / "taxi_zones.zip"47 shapes = read_zone_shapes(shapes_path) if shapes_path.exists() else {}48 xy, pos = stop_positions(inst, zones)49 top_ids = set(zones.LocationID)50 lid_of = dict(zip(zones.idx, zones.LocationID))51 fig, axes = plt.subplots(1, 2, figsize=(15, 8.2), gridspec_kw={"width_ratios": [1, 1.25]})52 has_jfk = any(zones.set_index("idx").loc[o.zone, "zone"] == "JFK Airport" for o in inst.orders)53 for ax, zoom in zip(axes, (False, True)):54 for lid, rings in shapes.items():55 hot = lid in top_ids56 for ring in rings:57 ax.fill(ring[:, 0], ring[:, 1], facecolor="#e9e8e4" if hot else "#f4f3f0",58 edgecolor="#c9c8c2" if hot else "#e6e5e1", linewidth=0.5, zorder=1)59 for v in inst.vans:60 r = res.routes.get(v.id)61 if not r or not r.sequence:62 continue63 pts = [xy[inst.depot]] + [pos[i] for i in r.sequence] + ([xy[inst.depot]] if inst.return_to_depot else [])64 for k, (a, b) in enumerate(zip(pts[:-1], pts[1:])):65 ax.annotate("", xy=b, xytext=a, zorder=3,66 arrowprops=dict(arrowstyle="-|>", color=P.VAN_COLORS[v.id], lw=1.8,67 shrinkA=6, shrinkB=6, alpha=0.9,68 linestyle="--" if k == len(pts) - 2 else "-"))69 for n, oid in enumerate(r.sequence, 1):70 late = r.lateness[n - 1] > 1e-971 ax.scatter(*pos[oid], s=150, color=P.VAN_COLORS[v.id], edgecolor=P.LATE if late else P.SURFACE,72 linewidth=2, zorder=4)73 ax.annotate(str(oid), pos[oid], ha="center", va="center", fontsize=7, color="white",74 weight="bold", zorder=5)75 ax.scatter(*xy[inst.depot], marker="*", s=420, color=P.INK, edgecolor=P.SURFACE, zorder=6)76 ax.set_aspect("equal"); ax.grid(False); ax.set_xticks([]); ax.set_yticks([])77 for s in ax.spines.values():78 s.set_visible(False)79 manh = xy[[i for i in range(len(xy)) if zones.set_index("idx").loc[i, "zone"] != "JFK Airport"]]80 if zoom:81 ax.set_xlim(manh[:, 0].min() - 1.0, manh[:, 0].max() + 1.0)82 ax.set_ylim(manh[:, 1].min() - 0.8, manh[:, 1].max() + 0.8)83 ax.set_title("Zoom: Manhattan (numbers = order ids; red ring = late)")84 if has_jfk:85 ax.annotate("legs to/from JFK Airport continue off-map", (0.98, 0.02), xycoords="axes fraction",86 ha="right", fontsize=8.5, color=P.INK2)87 for i in range(len(xy)):88 if zones.set_index("idx").loc[i, "zone"] == "JFK Airport":89 continue90 ax.annotate(zones.set_index("idx").loc[i, "zone"][:22], xy[i], xytext=(0, -14),91 textcoords="offset points", ha="center", fontsize=6.5, color=P.MUTED, zorder=2)92 else:93 allp = xy94 ax.set_xlim(allp[:, 0].min() - 3, allp[:, 0].max() + 3)95 ax.set_ylim(allp[:, 1].min() - 3, allp[:, 1].max() + 3)96 ax.set_title("Routes on the 20 busiest taxi zones" + (" (incl. JFK)" if has_jfk else ""))97 handles = [plt.Line2D([], [], color=c, lw=2.5, label=f"van {v}") for v, c in enumerate(P.VAN_COLORS[:len(inst.vans)])]98 handles += [plt.Line2D([], [], marker="*", color=P.INK, lw=0, ms=14, label="depot"),99 plt.Line2D([], [], color=P.INK2, lw=1.5, ls="--", label="return leg")]100 axes[0].legend(handles=handles, loc="lower left")101 return P.save(fig, path)102 103 104def draw_gantt(inst, res, path):105 rows = []106 for v in inst.vans:107 r = res.routes.get(v.id)108 if r:109 rows += [(v.id, k, oid) for k, oid in enumerate(r.sequence)]110 fig, ax = plt.subplots(figsize=(13, 0.36 * len(rows) + 1.8))111 y = 0112 yt, yl = [], []113 for v in inst.vans:114 r = res.routes.get(v.id)115 if not r:116 continue117 col = P.VAN_COLORS[v.id]118 ax.axhspan(y - 0.6, y + len(r.sequence) - 0.4, color=col, alpha=0.04, zorder=0)119 ax.annotate(f"van {v.id}\nshift {fmt_t(v.shift_start)}-{fmt_t(v.shift_end)}\nback {fmt_t(r.finish)}",120 (inst_min(inst) - 70, y + len(r.sequence) / 2 - 0.5), va="center", fontsize=8.5, color=P.INK2)121 for k, oid in enumerate(r.sequence):122 o = inst.order(oid)123 ax.barh(y, o.latest - o.earliest, left=o.earliest, height=0.62, color=col, alpha=0.22,124 edgecolor=col, linewidth=1, zorder=1)125 a, b = r.arrivals[k], r.begins[k]126 if b - a > 0.5:127 ax.plot([a, b], [y, y], color=P.INK2, lw=1, ls=":", zorder=2)128 late = r.lateness[k] > 1e-9129 ax.barh(y, o.service, left=b, height=0.36, color=P.LATE if late else col, zorder=3)130 ax.scatter([a], [y], marker="|", s=120, color=P.INK, zorder=4)131 if late:132 ax.annotate(f"late {r.lateness[k]:.0f} min", (b + o.service, y), xytext=(4, 0),133 textcoords="offset points", va="center", color=P.LATE, fontsize=8)134 yt.append(y); yl.append(f"#{oid} (zone {o.zone})")135 y += 1136 y += 0.6137 ax.set_yticks(yt); ax.set_yticklabels(yl, fontsize=8); ax.invert_yaxis()138 lo, hi = inst_min(inst) - 75, max(r.finish for r in res.routes.values() if r) + 15139 ticks = np.arange(60 * math.floor(lo / 60), hi + 1, 60)140 ax.set_xticks(ticks); ax.set_xticklabels([fmt_t(t) for t in ticks]); ax.set_xlim(lo, hi)141 ax.grid(axis="y", visible=False)142 ax.set_title("Timeline: time window (shaded), arrival (|), waiting (dotted), service (solid bar)")143 return P.save(fig, path)144 145 146def inst_min(inst):147 return min(v.shift_start for v in inst.vans)148 149 150def draw_constraint_graph(inst, res, path):151 csp = res.csp152 assign = res.assignment153 groups = {v.id: [i for i in sorted(assign) if assign[i] == v.id] for v in inst.vans}154 pos = {}155 nv = len(inst.vans)156 for v, ids in groups.items():157 c = 2 * np.pi * v / nv + np.pi / 2158 span = 2 * np.pi / nv * 0.78159 for k, i in enumerate(ids):160 t = c - span / 2 + span * (k + 0.5) / max(1, len(ids))161 pos[i] = np.array([np.cos(t), np.sin(t)])162 fig, ax = plt.subplots(figsize=(8, 8))163 for e, src in csp.edges.items():164 a, b = tuple(e)165 pa, pb = pos[a], pos[b]166 if src == "learned":167 ax.plot(*zip(pa, pb), color=P.LATE, lw=2.6, ls="--", zorder=2)168 mid = (pa + pb) / 2169 ax.annotate("learned", mid, color=P.LATE, fontsize=9, ha="center",170 bbox=dict(boxstyle="round,pad=0.2", fc=P.SURFACE, ec="none"))171 else:172 ax.plot(*zip(pa, pb), color=P.INK2, lw=1.4, zorder=1)173 for i, p in pos.items():174 ax.scatter(*p, s=520, color=P.VAN_COLORS[assign[i]], edgecolor=P.SURFACE, linewidth=2, zorder=3)175 ax.annotate(str(i), p, ha="center", va="center", color="white", weight="bold", zorder=4)176 for v in inst.vans:177 c = 2 * np.pi * v.id / nv + np.pi / 2178 ax.annotate(f"van {v.id}", 1.3 * np.array([np.cos(c), np.sin(c)]), ha="center", va="center",179 color=P.VAN_COLORS[v.id], weight="bold", fontsize=12)180 n_nn = sum(1 for s in csp.edges.values() if s == "nn")181 n_l = sum(1 for s in csp.edges.values() if s == "learned")182 handles = [plt.Line2D([], [], color=P.INK2, lw=1.4, label=f"different-van rule from NN times ({n_nn})"),183 plt.Line2D([], [], color=P.LATE, lw=2.6, ls="--", label=f"rule learned from a late route ({n_l})")]184 ax.legend(handles=handles, loc="lower center", bbox_to_anchor=(0.5, -0.08))185 ax.set_xlim(-1.5, 1.5); ax.set_ylim(-1.5, 1.5); ax.set_aspect("equal"); ax.axis("off")186 ax.set_title("CSP constraint graph: orders (nodes) grouped by assigned van")187 return P.save(fig, path)188 189 190def draw_stats(search_rows, csp_rows, loops, nn_rows, path):191 fig, axes = plt.subplots(2, 2, figsize=(13, 8.5))192 ax = axes[0, 0]193 algos = ["UCS", "Greedy", "A*"]194 vans = sorted({r["van"] for r in search_rows})195 for k, a in enumerate(algos):196 vals = [next(r["expanded"] for r in search_rows if r["van"] == v and r["algo"] == a) for v in vans]197 b = ax.bar(np.arange(len(vans)) + (k - 1) * 0.27, vals, 0.25, color=P.SERIES[k], label=a,198 edgecolor=P.SURFACE, linewidth=2)199 ax.bar_label(b, fontsize=8, color=P.INK2)200 ax.set_xticks(range(len(vans))); ax.set_xticklabels([f"van {v}" for v in vans])201 ax.set_title("Search: nodes expanded per van (this day)"); ax.legend()202 203 ax = axes[0, 1]204 labels = [r["solver"] for r in csp_rows]205 for k, key in enumerate(["assignments", "backtracks"]):206 b = ax.bar(np.arange(len(labels)) + (k - 0.5) * 0.38, [r[key] for r in csp_rows], 0.36,207 color=P.SERIES[k], label=key, edgecolor=P.SURFACE, linewidth=2)208 ax.bar_label(b, fontsize=8, color=P.INK2)209 ax.set_xticks(range(len(labels))); ax.set_xticklabels(labels, fontsize=8.5)210 ax.set_title("CSP pruning on a hard rush-hour instance (20 orders)"); ax.legend()211 212 ax = axes[1, 0]213 ls = [l["loop"] for l in loops]214 late = [l.get("late_stops", 0) for l in loops]215 learned = [l["n_edges_learned"] for l in loops]216 b = ax.bar(np.array(ls) - 0.2, late, 0.38, color=P.LATE, label="late routes after A*")217 ax.bar_label(b, fontsize=9, color=P.INK2)218 b = ax.bar(np.array(ls) + 0.2, learned, 0.38, color=P.SERIES[0], label="learned rules in CSP")219 ax.bar_label(b, fontsize=9, color=P.INK2)220 ax.set_xticks(ls); ax.set_xlabel("integration loop"); ax.set_title("Feedback loop: CSP <- late pair <- A*")221 ax.legend()222 223 ax = axes[1, 1]224 names = [r["model"] for r in nn_rows]225 vals = [r["late_mean"] for r in nn_rows]226 b = ax.bar(range(len(names)), vals, 0.6, color=P.SERIES[:len(names)], edgecolor=P.SURFACE, linewidth=2)227 ax.bar_label(b, labels=[f"{v:.2f}" for v in vals], fontsize=9, color=P.INK2)228 ax.set_xticks(range(len(names))); ax.set_xticklabels(names)229 ax.set_title("NN on vs off: late deliveries in replay of real test-week trips")230 ax.set_ylabel("late deliveries per day (mean of 50 replays)")231 return P.save(fig, path)232 233 234# ------------------------------------------------------------------------------ main235def main():236 ap = argparse.ArgumentParser()237 ap.add_argument("--seed", type=int, default=4)238 ap.add_argument("--impossible", action="store_true")239 ap.add_argument("--reps", type=int, default=50)240 args = ap.parse_args()241 C.OUTPUTS.mkdir(exist_ok=True)242 trips, zones, dist, mats = load_all()243 zname = zones.set_index("idx").zone244 inst = generate_instance(20, args.seed, n_zones=len(zones), zone_weights=zone_weights(zones),245 widths=(15, 30, 45, 60), depot=DEPOT_ZONE)246 if args.impossible:247 jfk = int(zones[zones.zone == "JFK Airport"].idx.iloc[0])248 o = inst.orders[0]249 o.zone, o.earliest, o.latest = jfk, 8 * 60, 8 * 60 + 20250 inst.name += "_impossible"251 tm = time_model(mats, "nn", inst.weekday)252 253 buf = io.StringIO()254 255 def out(*a):256 print(*a)257 print(*a, file=buf)258 259 out(f"=== INPUT: {inst.name}, weekday {C.WEEKDAY_NAMES[inst.weekday]}, depot zone {inst.depot} ({zname[inst.depot]})")260 out(f"travel times: NN-predicted taxi minutes x {C.VAN_FACTOR} (van factor)")261 for v in inst.vans:262 out(f" van {v.id}: shift {fmt_t(v.shift_start)}-{fmt_t(v.shift_end)}, capacity {v.capacity}, max {v.max_stops} stops")263 out(" order zone window demand service")264 for o in inst.orders:265 out(f" {o.id:5d} {o.zone:2d} {zname[o.zone][:28]:<28} {fmt_t(o.earliest)}-{fmt_t(o.latest)} {o.demand:5d} {o.service:5.0f}")266 267 traces = {}268 res = plan(inst, tm, traces=traces)269 out("\n=== INTEGRATION LOOP (NN -> CSP -> A*)")270 for l in res.loops:271 s = l["csp_stats"]272 out(f" loop {l['loop']}: CSP {l['csp_status']} (rules: {l['n_edges_nn']} from NN, {l['n_edges_learned']} learned; "273 f"{s['assignments']} assignments, {s['backtracks']} backtracks)")274 for nr in l["new_rules"]:275 out(f" van {nr['van']} late by {nr['late_by']:.1f} min -> new rule: orders {nr['pair'][0]} and {nr['pair'][1]} on different vans")276 if res.status == "no_solution":277 out("\n" + res.message)278 with open(C.OUTPUTS / "demo_impossible.txt", "w", encoding="utf-8") as f:279 f.write(buf.getvalue())280 return281 out("\n=== SOLUTION")282 out(describe_plan(inst, res))283 # loop-level late counts for the figure (re-plan each loop's assignment is costly; use recorded rules)284 for l in res.loops:285 l["late_stops"] = len(l["new_rules"])286 287 out("\n=== A* SEARCH TRACE (van 0, first 12 expansions; f = g + h, minutes)")288 for t in traces.get(0, [])[:12]:289 out(f" #{t['step']:3d} at {str(t['at']):>5} visited {str(t['visited']):<16} clock {fmt_t(t['clock'])} "290 f"g {t['g']:6.1f} h {t['h']:6.1f} f {t['f']:6.1f}")291 search_rows = []292 out("\n=== SEARCH COUNTERS per van (same stops, three algorithms)")293 for v in inst.vans:294 r = res.routes[v.id]295 if not r.sequence:296 continue297 prob = van_problem(inst, tm, v, r.sequence)298 for label, algo in (("UCS", "ucs"), ("Greedy", "greedy"), ("A*", "astar")):299 rr = solve_route(prob, algo)300 search_rows.append({"van": v.id, "algo": label, "expanded": rr.stats["expanded"],301 "generated": rr.stats["generated"], "cost": rr.cost})302 out(f" van {v.id} {label:<6} expanded {rr.stats['expanded']:5d} generated {rr.stats['generated']:5d} "303 f"route time {rr.cost:6.1f} min {1000*rr.stats['time_s']:.1f} ms")304 305 out("\n=== CSP PRUNING: same rules, solvers switched on/off")306 out(" (a) this day's constraint set")307 for fc, ac3, label in ((False, False, "BT"), (True, False, "BT+FC"), (True, True, "BT+FC+AC3")):308 s = res.csp.solve(fc=fc, ac3=ac3, var_order="input").stats309 out(f" {label:<10} assignments {s['assignments']:5d} backtracks {s['backtracks']:5d} checks {s['checks']:5d}")310 hard = generate_instance(20, 4, n_zones=len(zones), zone_weights=zone_weights(zones), cap_slack=1.0,311 max_stops=8, day=(540, 720), widths=(15, 30, 45))312 hcsp = AssignmentCSP(hard, time_model(mats, "nn", hard.weekday))313 out(f" (b) hard rush-hour instance: 20 orders, windows 09:00-12:00, zero capacity slack, {len(hcsp.edges)} NN rules")314 csp_rows = []315 for fc, ac3, label, vo in ((False, False, "BT", "input"), (True, False, "BT+FC", "input"),316 (True, True, "BT+FC+AC3", "input"), (True, True, "+degree order", "degree")):317 s = hcsp.solve(fc=fc, ac3=ac3, var_order=vo).stats318 csp_rows.append({**s, "solver": label})319 out(f" {label:<14} assignments {s['assignments']:5d} backtracks {s['backtracks']:5d} "320 f"fc_prunes {s['fc_prunes']:5d} ac3_prunes {s['ac3_prunes']:4d} {1000*s['time_s']:.1f} ms")321 322 out(f"\n=== NN ON vs OFF: replay with real held-out trip durations ({args.reps} replays, test week)")323 sampler = ReplaySampler(trips[trips.split == "test"])324 nn_rows = []325 for kind in ("nn", "hist", "dist"):326 tmk = time_model(mats, kind, inst.weekday)327 rk = res if kind == "nn" else plan(inst, tmk, suggest=False)328 rp = replay(inst, rk, sampler, tmk, n_reps=args.reps, seed=0)329 label = {"nn": "NN on", "hist": "hist-avg", "dist": "NN off (distance)"}[kind]330 nn_rows.append({"model": label, **rp})331 out(f" {label:<18} plan {rk.status:<4} late per day {rp['late_mean']:.2f} ± {rp['late_std']:.2f} "332 f"({100*rp['late_frac']:.1f}%), lateness {rp['late_minutes_mean']:.1f} min")333 334 figs = [draw_map(inst, res, zones, C.OUTPUTS / "demo_map.png"),335 draw_gantt(inst, res, C.OUTPUTS / "demo_gantt.png"),336 draw_constraint_graph(inst, res, C.OUTPUTS / "demo_constraint_graph.png"),337 draw_stats(search_rows, csp_rows, res.loops, nn_rows, C.OUTPUTS / "demo_stats.png")]338 out("\nfigures: " + ", ".join(str(f.relative_to(C.ROOT)) for f in figs))339 write_html(buf.getvalue(), figs)340 out("report: outputs/demo_report.html")341 342 343def write_html(text, figs):344 imgs = "\n".join(f'<figure><img src="{f.name}" alt="{f.stem}"><figcaption>{f.stem}</figcaption></figure>' for f in figs)345 extra = [p for p in ("time_matrix_nn.png", "nn_mae_by_hour.png", "exp1_search.png", "exp2_csp.png",346 "exp3_lateness.png") if (C.OUTPUTS / p).exists()]347 imgs += "\n".join(f'<figure><img src="{p}" alt="{p}"><figcaption>{p}</figcaption></figure>' for p in extra)348 page = f"""<!doctype html><html><head><meta charset="utf-8"><title>Delivery Routing Demo</title>349<meta name="viewport" content="width=device-width, initial-scale=1">350<style>351:root{{--bg:#fcfcfb;--ink:#0b0b0b;--ink2:#52514e;--line:#e4e3df}}352body{{background:var(--bg);color:var(--ink);font:15px/1.5 system-ui,sans-serif;max-width:1200px;margin:0 auto;padding:16px}}353h1{{font-size:24px}} pre{{background:#f4f3f0;border:1px solid var(--line);padding:12px;overflow-x:auto;font-size:12px}}354img{{max-width:100%;border:1px solid var(--line)}} figcaption{{color:var(--ink2);font-size:13px}} figure{{margin:24px 0}}355</style></head><body>356<h1>P4 Delivery Routing with Learned Travel Times: demo run</h1>357<p>NN travel-time matrix → CSP assigns orders to vans → A* orders each van's stops → late pair becomes a new CSP rule.</p>358{imgs}359<h2>Console output</h2><pre>{html.escape(text)}</pre></body></html>"""360 (C.OUTPUTS / "demo_report.html").write_text(page, encoding="utf-8")361 362 363if __name__ == "__main__":364 main()