#!/usr/bin/env python3 """No.04 Braess のパラドックス — 車を1台ずつ乗り換えさせる実験器（E1〜E4）。 網は Easley & Kleinberg (2010) 第8章 図8.1・8.2 の例をそのまま使う: A→C: x/100 分（x はその区間を通る台数） C→B: 45 分 A→D: 45 分 D→B: x/100 分 足す道 C→D: 0 分（台数によらない） 経路は 上回り A→C→B ／ 下回り A→D→B ／ 近道 A→C→D→B（近道は道を足した世界だけ）。 動き（全実験共通）: 最初、各車は上回りか下回りを五分五分で選ぶ（乱数）。その後、毎回 でたらめに選んだ1台が、いまの混み方で一番早く着く経路を計算し、今より早くなるときだけ 乗り換える（同じ時間なら乗り換えない）。誰も乗り換えたくない状態（均衡）になったら止まる。 「元の世界」（道なし）と「道を足した世界」は、同じ乱数の種で回す＝最初の分かれ方も、 選ばれる車の順番も同じ。違うのは C→D の道が在るかどうかだけ。 時間は整数（1/100 分＝センチ分）で計算する。x/100 分 = x センチ分、45 分 = 4500 センチ分。 標準ライブラリのみ。予想（predictions.json）は本番の種を回す前に凍結すること。 """ from __future__ import annotations import argparse import datetime as dt import json import platform import random import statistics import time from pathlib import Path from typing import Any UPPER, LOWER, SHORT = 0, 1, 2 ROUTE_NAMES = ("upper", "lower", "short") FIXED = 4500 # 45 分（センチ分） DEFAULT_CARS = 4000 DEFAULT_SEEDS = list(range(4001, 5001)) SWEEP_CARS = [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000] SWEEP_SEEDS = 100 # 掃引は本番の種の先頭 100 個 MAX_STEPS = 10_000_000 EARLY_SNAPS = 40 # 道を足した世界の写し: 最初の 40 回の乗り換えは毎回 LATE_EVERY = 40 # その後は 40 回ごと def route_times(u: int, l: int, m: int) -> tuple[int, int, int]: """上回り・下回り・近道の所要時間（センチ分）。s=A→C の台数、t=D→B の台数。""" s, t = u + m, l + m return s + FIXED, FIXED + t, s + t def total_time(u: int, l: int, m: int) -> int: tu, tl, tm = route_times(u, l, m) return u * tu + l * tl + m * tm def time_after_move(src: int, dst: int, counts: list[int]) -> int: """src の車が1台 dst へ移ったときの、その車の所要時間（センチ分）。""" c = counts[:] c[src] -= 1 c[dst] += 1 return route_times(*c)[dst] def best_move(src: int, counts: list[int], road: bool) -> int | None: """src の車にとって今より早い経路（無ければ None）。同じ時間の候補は 上・下・近道 の順。""" current = route_times(*counts)[src] best, best_t = None, current for dst in (UPPER, LOWER, SHORT): if dst == src or (dst == SHORT and not road): continue t = time_after_move(src, dst, counts) if t < best_t: best, best_t = dst, t return best def is_equilibrium(counts: list[int], road: bool) -> bool: return all(counts[r] == 0 or best_move(r, counts, road) is None for r in (UPPER, LOWER, SHORT)) def run_world(seed: int, road: bool, cars: int = DEFAULT_CARS, snaps: str | None = None, max_steps: int = MAX_STEPS, open_after_eq: bool = False) -> dict[str, Any]: """1つの世界を均衡まで回す。snaps='all' なら乗り換えごと、'road' なら早い 40 回＋40 回ごとの写し。 open_after_eq=True（見せる回だけ・2026-11-05 追加）: まず道なしで均衡まで回し、同じ乱数の続きで C→D の道を開けて、そこから均衡まで回す（「道なしの世界が落ち着いた後に、道を足した」場合）。 steps・moves・写し・最初の乗り換えなどは、道を開けた後だけを数える。E1 の組（同じ最初の分かれ方から 両方の世界を回す）とは別の見せ方で、予想（predictions.json）の対象ではない。 """ rng = random.Random(seed) route = [UPPER if rng.random() < 0.5 else LOWER for _ in range(cars)] counts = [route.count(UPPER), route.count(LOWER), 0] init = counts[:] pre_steps, opened_at = 0, None if open_after_eq: if not road: raise ValueError("open_after_eq は道を足す世界だけ") while not is_equilibrium(counts, False): if pre_steps >= max_steps: raise RuntimeError("道なしの段が均衡に着かない") pre_steps += 1 i = rng.randrange(cars) src = route[i] dst = best_move(src, counts, False) if dst is None: continue route[i] = dst counts[src] -= 1 counts[dst] += 1 opened_at = counts[:] steps = moves = 0 first_switch = None # 最初に近道へ移った車の、移った直後の所要時間（センチ分） tot0 = total_time(*counts) min_total, min_total_at = tot0, 0 # 乗り換えの途中で、全体の合計時間が一番短かった時（moves の値） snap_list: list[dict[str, Any]] = [] def save() -> None: tu, tl, tm = route_times(*counts) snap_list.append({"step": steps, "moves": moves, "u": counts[0], "l": counts[1], "m": counts[2], "t_upper": tu, "t_lower": tl, "t_short": tm, "total": total_time(*counts)}) if snaps: save() eq = is_equilibrium(counts, road) while not eq: if steps >= max_steps: break steps += 1 i = rng.randrange(cars) src = route[i] dst = best_move(src, counts, road) if dst is None: continue route[i] = dst counts[src] -= 1 counts[dst] += 1 moves += 1 if dst == SHORT and first_switch is None: first_switch = route_times(*counts)[SHORT] tot = total_time(*counts) if tot < min_total: min_total, min_total_at = tot, moves eq = is_equilibrium(counts, road) if snaps == "all" or (snaps == "road" and (moves <= EARLY_SNAPS or moves % LATE_EVERY == 0 or eq)): save() tu, tl, tm = route_times(*counts) used = [r for r in (UPPER, LOWER, SHORT) if counts[r] > 0] times_used = sorted({(tu, tl, tm)[r] for r in used}) out = { "seed": seed, "road": road, "cars": cars, "init": init, "final": counts[:], "steps": steps, "moves": moves, "equilibrium": eq, "t_upper": tu, "t_lower": tl, "t_short": tm if road else None, "times_used": times_used, "total": total_time(*counts), "first_switch": first_switch, "min_total": min_total, "min_total_at": min_total_at, # 近道の車が1台だけ上回りへ戻ったときの所要時間（道を足した世界だけ） "back_to_upper": time_after_move(SHORT, UPPER, counts) if road and counts[SHORT] > 0 else None, } if open_after_eq: out.update({"open_after_eq": True, "pre_steps": pre_steps, "opened_at": opened_at}) if snaps: out["snaps"] = snap_list return out def median_low(xs: list[int]) -> int: return statistics.median_low(xs) def p90(xs: list[int]) -> int: ys = sorted(xs) k = -(-9 * len(ys) // 10) # ceil(0.9 n) return ys[k - 1] def theory_minutes(cars: int, road: bool) -> float: """連続の流れ（Wardrop 均衡）の閉じた形。式は Easley & Kleinberg (2010) 第8章の考え方を台数で一般化。""" if not road: return cars / 200 + 45 if cars <= 4500: # 全員が近道 return cars / 50 if cars <= 9000: # 3経路とも使い、どれも 90 分 return 90.0 return cars / 200 + 45 # 近道は使われない def social_optimum(cars: int) -> dict[str, Any]: """道を足した世界で、全員の合計時間が一番短くなる分け方（整数・全部の m を調べ、上下は半々に近い2通り）。""" best = None for m in range(cars + 1): rest = cars - m for u in {rest // 2, rest - rest // 2}: l = rest - u tot = total_time(u, l, m) if best is None or tot < best[0]: best = (tot, u, l, m) tot, u, l, m = best return {"total": tot, "avg_minutes": tot / cars / 100, "u": u, "l": l, "m": m, "method": "m を 0〜台数まで全部、上下はその残りを半々に近い2通りで調べた（合計は u と l の差について下に凸なので半々が最小）"} def summarize_pair(base: list[dict[str, Any]], road: list[dict[str, Any]], cars: int) -> dict[str, Any]: n = len(base) bm = [r["total"] for r in base] rm = [r["total"] for r in road] base_min = [r["total"] / cars / 100 for r in base] road_min = [r["total"] / cars / 100 for r in road] worse = sum(1 for b, r in zip(bm, rm) if r > b) better = sum(1 for b, r in zip(bm, rm) if r < b) road_steps = [r["steps"] for r in road] base_steps = [r["steps"] for r in base] return { "n": n, "base_equilibrium": sum(r["equilibrium"] for r in base), "road_equilibrium": sum(r["equilibrium"] for r in road), "base_minutes_set": sorted(set(base_min)), "road_minutes_set": sorted(set(road_min)), "base_minutes": base_min[0] if len(set(base_min)) == 1 else None, "road_minutes": road_min[0] if len(set(road_min)) == 1 else None, "increase_minutes": (road_min[0] - base_min[0]) if len(set(base_min)) == 1 and len(set(road_min)) == 1 else None, "worse": worse, "better": better, "same": n - worse - better, "road_all_short": sum(1 for r in road if r["final"][SHORT] == cars), "base_balanced": sum(1 for r in base if r["final"][UPPER] == r["final"][LOWER]), "road_moves_set": sorted({r["moves"] for r in road}), "road_moves_all_cars": sum(1 for r in road if r["moves"] == cars), "base_moves_median": median_low([r["moves"] for r in base]), "base_moves_max": max(r["moves"] for r in base), "base_steps_median": median_low(base_steps), "base_steps_max": max(base_steps), "road_steps_median": median_low(road_steps), "road_steps_p90": p90(road_steps), "road_steps_min": min(road_steps), "road_steps_max": max(road_steps), "road_steps_mean": sum(road_steps) / n, "first_switch_minutes_min": min(r["first_switch"] for r in road) / 100, "first_switch_minutes_max": max(r["first_switch"] for r in road) / 100, "road_min_avg_minutes_min": min(r["min_total"] for r in road) / cars / 100, "road_min_avg_minutes_max": max(r["min_total"] for r in road) / cars / 100, "road_min_below_base": sum(1 for b, r in zip(base, road) if r["min_total"] < b["total"]), "back_to_upper_minutes_set": sorted({r["back_to_upper"] / 100 for r in road}), "init_upper_min": min(r["init"][UPPER] for r in base), "init_upper_max": max(r["init"][UPPER] for r in base), "steps": {"base": base_steps, "road": road_steps}, } def sweep(seeds: list[int], cars_list: list[int]) -> dict[str, Any]: by = {} for cars in cars_list: rows = {} for road in (False, True): runs = [run_world(s, road, cars) for s in seeds] mins = [r["total"] / cars / 100 for r in runs] rows["road" if road else "base"] = { "n": len(runs), "equilibrium": sum(r["equilibrium"] for r in runs), "minutes_mean": sum(mins) / len(mins), "minutes_min": min(mins), "minutes_max": max(mins), "short_mean": sum(r["final"][SHORT] for r in runs) / len(runs), "theory": theory_minutes(cars, road), } b, r = rows["base"], rows["road"] by[str(cars)] = {"cars": cars, "base": b, "road": r, "diff_mean": r["minutes_mean"] - b["minutes_mean"], "road_worse": r["minutes_min"] > b["minutes_max"] + 1e-9, "road_better": r["minutes_max"] < b["minutes_min"] - 1e-9} devs = [abs(x[w]["minutes_mean"] - x[w]["theory"]) for x in by.values() for w in ("base", "road")] worse = [x["cars"] for x in by.values() if x["road_worse"]] better = [x["cars"] for x in by.values() if x["road_better"]] return {"cars": cars_list, "seeds_first": seeds[0], "seeds_last": seeds[-1], "n_per_point": len(seeds), "by_cars": by, "max_abs_dev": max(devs), "worse_cars": worse, "better_cars": better, "same_cars": [c for c in cars_list if c not in worse and c not in better]} def build_summary(seeds: list[int], cars: int = DEFAULT_CARS, sweep_n: int = SWEEP_SEEDS, sweep_cars: list[int] | None = None) -> dict[str, Any]: t0 = time.monotonic() base = [run_world(s, False, cars) for s in seeds] road = [run_world(s, True, cars) for s in seeds] pair = summarize_pair(base, road, cars) fs = seeds[0] fb = run_world(fs, False, cars, snaps="all") fr = run_world(fs, True, cars, snaps="road") fo = run_world(fs, True, cars, snaps="road", open_after_eq=True) sw = sweep(seeds[:sweep_n], sweep_cars or SWEEP_CARS) opt = social_optimum(cars) theory = { "before_minutes": theory_minutes(cars, False), "after_minutes": theory_minutes(cars, True), "increase_minutes": theory_minutes(cars, True) - theory_minutes(cars, False), "back_to_upper_minutes": cars / 100 + 45, "first_switch_minutes_balanced": route_times(cars // 2 - 1, cars // 2, 1)[SHORT] / 100, "worse_range": [3000, 9000], "optimum": opt, "textbook": {"before": 65, "after": 80, "source": "Easley & Kleinberg (2010) 第8章 図8.1・8.2（4,000台）"}, "method": "閉じた形（連続の流れの均衡）＝道なし 台数/200+45・道あり 台数≦4500 で 台数/50、4500〜9000 で 90、9000 以上で 台数/200+45", } return { "model": "braess_best_response", "config": {"cars": cars, "seeds": seeds, "n_per_condition": len(seeds), "fixed_minutes": FIXED // 100, "variable": "x/100", "variable_divisor": 100, "new_road_minutes": 0, "unit": "centi-minute", "early_snaps": EARLY_SNAPS, "late_every": LATE_EVERY, "max_steps": MAX_STEPS, "sweep_cars": sweep_cars or SWEEP_CARS, "sweep_n": sweep_n, "python": platform.python_version(), "rng": "random.Random(seed)"}, "groups": {"pair": {"base": {"n": len(base)}, "road": {"n": len(road)}}, "sweep": {str(c): {"n": sweep_n} for c in (sweep_cars or SWEEP_CARS)}}, "pair": pair, "featured": {"seed": fs, "rule": "本番の種の最初の1つ", "base": fb, "road": fr, "road_open": fo, "road_open_rule": "同じ種で、道なしの世界が均衡で止まった後に、乱数の続きで C→D を開けた回（ショートの見せ方・予想の対象外）"}, "sweep": sw, "theory": theory, "generated_at": dt.datetime.now(dt.timezone(dt.timedelta(hours=9))).isoformat(timespec="seconds"), "runtime_s": round(time.monotonic() - t0, 3), } def metric_value(summary: dict[str, Any], metric: str) -> Any: value: Any = summary for part in metric.split("."): value = value[int(part)] if isinstance(value, list) else value[part] return value def prediction_hit(prediction: dict[str, Any], observed: Any) -> bool: if observed is None: return False kind = prediction["kind"] if kind == "eq": return observed == prediction["value"] if kind == "bool": return isinstance(observed, bool) and observed is prediction["value"] if kind == "range": return prediction["lo"] <= observed <= prediction["hi"] if kind == "max": return observed <= prediction["hi"] if kind == "min": return observed >= prediction["lo"] if kind == "in": return observed in prediction["values"] raise ValueError(f"unknown prediction kind: {kind}") def check_predictions(summary: dict[str, Any], predictions: dict[str, Any] | Path) -> dict[str, Any]: if isinstance(predictions, Path): predictions = json.loads(predictions.read_text(encoding="utf-8")) checks = [] for p in predictions["predictions"]: observed = metric_value(summary, p["metric"]) text = f"{observed:.10g}" if isinstance(observed, float) else str(observed) checks.append({"id": p["id"], "text": p["text"], "metric": p["metric"], "observed": observed, "observed_text": text, "hit": prediction_hit(p, observed)}) return {"n": len(checks), "hits": sum(c["hit"] for c in checks), "checks": checks} def report_md(s: dict[str, Any], pc: dict[str, Any] | None) -> str: p, sw, th, c = s["pair"], s["sweep"], s["theory"], s["config"] rows = [ "# ひとつ変えてみる No.04 — Braess のパラドックス 結果表", "", f"網: Easley & Kleinberg (2010) 第8章 図8.1・8.2（{c['cars']:,}台・A→C と D→B は 台数/100 分・C→B と A→D は 45 分・足す道 C→D は 0 分）。", f"動き: でたらめに選んだ1台が、今より早くなる経路へ乗り換える。誰も乗り換えたくなくなったら止まる。種 {c['seeds'][0]}〜{c['seeds'][-1]}（n={p['n']}）。", "", "| 項目 | 道なし（元の世界） | 道を足した世界 |", "|---|---:|---:|", f"| 均衡で止まった回 | {p['base_equilibrium']}/{p['n']} | {p['road_equilibrium']}/{p['n']} |", f"| 止まった時の所要時間（分） | {', '.join(f'{x:.2f}' for x in p['base_minutes_set'])} | {', '.join(f'{x:.2f}' for x in p['road_minutes_set'])} |", f"| 全員が近道を使った回 | — | {p['road_all_short']}/{p['n']} |", f"| 乗り換えの回数 | 中央値 {p['base_moves_median']}・最大 {p['base_moves_max']} | {', '.join(str(x) for x in p['road_moves_set'])}（全回） |", f"| 止まるまでに選ばれた回数 | 中央値 {p['base_steps_median']}・最大 {p['base_steps_max']} | 中央値 {p['road_steps_median']:,}・p90 {p['road_steps_p90']:,}・最大 {p['road_steps_max']:,} |", "", f"- 道を足して遅くなった回 {p['worse']}・早くなった回 {p['better']}・同じ {p['same']}（n={p['n']}）。", f"- 最初に近道へ移った車の所要時間: {p['first_switch_minutes_min']:.2f}〜{p['first_switch_minutes_max']:.2f} 分。", f"- 途中で全体の平均が一番短かった時: {p['road_min_avg_minutes_min']:.4f}〜{p['road_min_avg_minutes_max']:.4f} 分（道なしの止まった値より短くなった回 {p['road_min_below_base']}/{p['n']}）。", f"- 止まった後、近道の車が1台だけ上回りへ戻ると: {', '.join(f'{x:.2f}' for x in p['back_to_upper_minutes_set'])} 分。", f"- 理論（閉じた形）: 道なし {th['before_minutes']:.2f}・道あり {th['after_minutes']:.2f}・戻ると {th['back_to_upper_minutes']:.2f} 分。教科書の値 {th['textbook']['before']}→{th['textbook']['after']} 分。", f"- みんなで相談して分けた場合（道あり・合計が最小）: 平均 {th['optimum']['avg_minutes']:.4f} 分（上 {th['optimum']['u']:,}・下 {th['optimum']['l']:,}・近道 {th['optimum']['m']:,}）。", "", f"## 台数を変えた掃引（各 n={sw['n_per_point']}・種 {sw['seeds_first']}〜{sw['seeds_last']}）", "", "| 台数 | 道なし（分） | 道あり（分） | 理論 道なし | 理論 道あり | 近道の台数（平均） |", "|---:|---:|---:|---:|---:|---:|", ] for x in sw["by_cars"].values(): b, r = x["base"], x["road"] rows.append(f"| {x['cars']:,} | {b['minutes_mean']:.3f} | {r['minutes_mean']:.3f} | {b['theory']:.2f} | {r['theory']:.2f} | {r['short_mean']:.1f} |") rows += ["", f"- 道を足して遅くなった台数: {', '.join(f'{c:,}' for c in sw['worse_cars'])}／早くなった台数: {', '.join(f'{c:,}' for c in sw['better_cars'])}／同じ: {', '.join(f'{c:,}' for c in sw['same_cars'])}。", f"- 理論とのずれの最大: {sw['max_abs_dev']:.4f} 分。"] if pc: rows += ["", f"## 予想（事前登録 {pc['n']}個）: 当たり {pc['hits']}・外れ {pc['n'] - pc['hits']}", "", "| id | 予想 | 実測 | 判定 |", "|---|---|---|---|"] rows += [f"| {x['id']} | {x['text']} | {x['observed_text']} | {'当たり' if x['hit'] else '外れ'} |" for x in pc["checks"]] return "\n".join(rows) + "\n" def dump_json(path: Path, value: dict[str, Any]) -> None: path.parent.mkdir(parents=True, exist_ok=True) path.write_text(json.dumps(value, ensure_ascii=False, indent=1) + "\n", encoding="utf-8") def main(argv: list[str] | None = None) -> int: ap = argparse.ArgumentParser() ap.add_argument("--results-dir", type=Path, help="summary.json・predictions_check.json の置き場") ap.add_argument("--report", type=Path, help="結果表（Markdown）の出力先") ap.add_argument("--seeds", type=int, nargs="+", default=DEFAULT_SEEDS) ap.add_argument("--cars", type=int, default=DEFAULT_CARS) ap.add_argument("--sweep-n", type=int, default=SWEEP_SEEDS) ap.add_argument("--no-predictions", action="store_true", help="予想の判定をしない（試し回し用）") a = ap.parse_args(argv) s = build_summary(a.seeds, a.cars, a.sweep_n) pc = None if not a.no_predictions: pc = check_predictions(s, Path(__file__).with_name("predictions.json")) if a.results_dir: dump_json(a.results_dir / "summary.json", s) if pc: dump_json(a.results_dir / "predictions_check.json", pc) if a.report: a.report.parent.mkdir(parents=True, exist_ok=True) a.report.write_text(report_md(s, pc), encoding="utf-8") p = s["pair"] print(json.dumps({"n": p["n"], "base_minutes": p["base_minutes"], "road_minutes": p["road_minutes"], "worse": p["worse"], "road_steps_median": p["road_steps_median"], "sweep_worse": s["sweep"]["worse_cars"], "hits": pc and pc["hits"], "runtime_s": s["runtime_s"]}, ensure_ascii=False)) return 0 if __name__ == "__main__": raise SystemExit(main())