diff --git a/tests/test_cpp_mode.py b/tests/test_cpp_mode.py index 9fde209..13ee817 100644 --- a/tests/test_cpp_mode.py +++ b/tests/test_cpp_mode.py @@ -5485,7 +5485,7 @@ def create_World(cpp=True): ################################# # DUE - solver_DUE = DTAsolvers.SolverDUE(create_World) + solver_DUE = DTAsolvers.SolverDUE(create_World, cpp=True) solver_DUE.solve(max_iter=20) # max_iter should be larger (e.g., 100). this is just for demonstration W_DUE = solver_DUE.W_sol W_DUE.analyzer.print_simple_stats(force_print=True) @@ -5509,7 +5509,7 @@ def create_World(cpp=True): ################################# # DSO by day-to-day: this is recommended - solver_DSO_D2D= DTAsolvers.SolverDSO_D2D(create_World) + solver_DSO_D2D= DTAsolvers.SolverDSO_D2D(create_World, cpp=True) solver_DSO_D2D.solve(max_iter=20) # max_iter should be larger (e.g., 100). this is just for demonstration W_DSO_D2D = solver_DSO_D2D.W_sol W_DSO_D2D.analyzer.print_simple_stats(force_print=True) @@ -5598,7 +5598,7 @@ def create_World(cpp=True): ################################# # DUE - solver_DUE = SolverDUE(create_World) + solver_DUE = SolverDUE(create_World, cpp=True) solver_DUE.solve(max_iter=20) # max_iter should be larger (e.g., 100). this is just for demonstration W_DUE = solver_DUE.W_sol W_DUE.analyzer.print_simple_stats(force_print=True) @@ -5706,7 +5706,7 @@ def create_World(cpp=True): ] } - solver_DUE = DTAsolvers.SolverDUE(create_World) + solver_DUE = DTAsolvers.SolverDUE(create_World, cpp=True) solver_DUE.solve(max_iter=20, route_sets=route_sets) W_DUE = solver_DUE.W_sol W_DUE.analyzer.print_simple_stats(force_print=True) @@ -5717,7 +5717,7 @@ def create_World(cpp=True): assert df_DUE_link["traffic_volume"][df_DUE_link["link"]=="offramp"].item() == 0 - solver_DSO_D2D= DTAsolvers.SolverDSO_D2D(create_World) + solver_DSO_D2D= DTAsolvers.SolverDSO_D2D(create_World, cpp=True) solver_DSO_D2D.solve(max_iter=20, route_sets=route_sets) W_DSO_D2D = solver_DSO_D2D.W_sol W_DSO_D2D.analyzer.print_simple_stats(force_print=True) @@ -5775,7 +5775,7 @@ def toll(t): return W # DUE - solver_DUE = SolverDUE(create_World) + solver_DUE = SolverDUE(create_World, cpp=True) solver_DUE.solve(max_iter=40, print_progress=False) W_DUE = solver_DUE.W_sol W_DUE.analyzer.print_simple_stats(force_print=False) @@ -7038,7 +7038,7 @@ def create_World(): W.adddemand("A", "D", 0, 2000, 0.4) return W - solver = SolverDUE(create_World) + solver = SolverDUE(create_World, cpp=True) solver.solve(max_iter=3, print_progress=False) W_sol = solver.W_sol W_sol.analyzer.print_simple_stats(force_print=False) @@ -7102,7 +7102,7 @@ def toll(t): W.adddemand("1", "3", 0, 3000, 0.3) return W - solver = SolverDUE(create_World) + solver = SolverDUE(create_World, cpp=True) solver.solve(max_iter=3, print_progress=False) W_sol = solver.W_sol assert W_sol.analyzer.total_travel_time > 0 @@ -7932,3 +7932,226 @@ def test_coverage_print_mode_enabled(): W.adddemand("A", "B", 0, 100, 0.3) W.exec_simulation() assert len(W.VEHICLES) > 0 + + +# ====================================================================== +# DTA solver tests for C++ mode (SolverDUE, SolverDSO_D2D) +# Based on 9x9 grid scenario from devlog/benchtestdta.py +# ====================================================================== + +def _create_dta_grid_world(seed, cpp): + """Helper: build a 9x9 grid DTA scenario (from benchtestdta.py).""" + W = uxsim.World( + name="", deltan=10, tmax=4800, + duo_update_time=300, + print_mode=0, save_mode=0, show_mode=0, + random_seed=seed, cpp=cpp, + ) + imax, jmax = 9, 9 + id_center = 4 + nodes = {} + for i in range(imax): + for j in range(jmax): + nodes[i, j] = W.addNode(f"n{(i, j)}", i, j, flow_capacity=1.6) + for i in range(imax): + for j in range(jmax): + free_flow_speed = 10 + if i != imax - 1: + if j == id_center: + free_flow_speed = 20 + W.addLink(f"l{(i, j, i+1, j)}", nodes[i, j], nodes[i+1, j], + length=1000, free_flow_speed=free_flow_speed) + free_flow_speed = 10 + if i != 0: + if j == id_center: + free_flow_speed = 20 + W.addLink(f"l{(i, j, i-1, j)}", nodes[i, j], nodes[i-1, j], + length=1000, free_flow_speed=free_flow_speed) + free_flow_speed = 10 + if j != jmax - 1: + if i == id_center: + free_flow_speed = 20 + W.addLink(f"l{(i, j, i, j+1)}", nodes[i, j], nodes[i, j+1], + length=1000, free_flow_speed=free_flow_speed) + free_flow_speed = 10 + if j != 0: + if i == id_center: + free_flow_speed = 20 + W.addLink(f"l{(i, j, i, j-1)}", nodes[i, j], nodes[i, j-1], + length=1000, free_flow_speed=free_flow_speed) + + demand_flow = 0.08 + demand_duration = 2400 + outer_ids = 3 + for n1 in [(0, j) for j in range(outer_ids, jmax - outer_ids)]: + for n2 in [(imax - 1, j) for j in range(outer_ids, jmax - outer_ids)]: + W.adddemand(nodes[n1], nodes[n2], 0, demand_duration, demand_flow) + for n2 in [(i, jmax - 1) for i in range(outer_ids, imax - outer_ids)]: + W.adddemand(nodes[n1], nodes[n2], 0, demand_duration, demand_flow) + for n1 in [(i, 0) for i in range(outer_ids, imax - outer_ids)]: + for n2 in [(i, jmax - 1) for i in range(outer_ids, imax - outer_ids)]: + W.adddemand(nodes[n1], nodes[n2], 0, demand_duration, demand_flow) + for n2 in [(imax - 1, j) for j in range(outer_ids, jmax - outer_ids)]: + W.adddemand(nodes[n1], nodes[n2], 0, demand_duration, demand_flow) + + return W + + +def test_dta_solver_due_grid_cpp(): + """SolverDUE runs on a 9x9 grid with cpp=True and returns positive TTT.""" + from uxsim.DTAsolvers import SolverDUE + + def create_World(): + return _create_dta_grid_world(seed=42, cpp=True) + + solver = SolverDUE(create_World, cpp=True) + solver.solve(max_iter=1, print_progress=False) + W_sol = solver.W_sol + W_sol.analyzer.print_simple_stats(force_print=False) + assert W_sol.analyzer.total_travel_time > 0 + assert W_sol.analyzer.trip_completed > 0 + + +def test_dta_solver_dso_d2d_grid_cpp(): + """SolverDSO_D2D runs on a 9x9 grid with cpp=True and returns positive TTT.""" + from uxsim.DTAsolvers import SolverDSO_D2D + + def create_World(): + return _create_dta_grid_world(seed=42, cpp=True) + + solver = SolverDSO_D2D(create_World, cpp=True) + solver.solve(max_iter=1, print_progress=False) + W_sol = solver.W_sol + W_sol.analyzer.print_simple_stats(force_print=False) + assert W_sol.analyzer.total_travel_time > 0 + assert W_sol.analyzer.trip_completed > 0 + + +@pytest.mark.flaky(reruns=5) +def test_dta_solver_due_grid_cpp_vs_python(): + """SolverDUE on 9x9 grid: C++ and Python produce similar TTT.""" + from uxsim.DTAsolvers import SolverDUE + + solver_cpp = SolverDUE(lambda: _create_dta_grid_world(seed=42, cpp=True), cpp=True) + solver_cpp.solve(max_iter=1, print_progress=False) + ttt_cpp = solver_cpp.W_sol.analyzer.total_travel_time + trips_cpp = solver_cpp.W_sol.analyzer.trip_completed + + solver_py = SolverDUE(lambda: _create_dta_grid_world(seed=42, cpp=False)) + solver_py.solve(max_iter=1, print_progress=False) + ttt_py = solver_py.W_sol.analyzer.total_travel_time + trips_py = solver_py.W_sol.analyzer.trip_completed + + assert trips_cpp == trips_py + assert ttt_cpp == pytest.approx(ttt_py, rel=0.2) + + +@pytest.mark.flaky(reruns=5) +def test_dta_solver_dso_d2d_grid_cpp_vs_python(): + """SolverDSO_D2D on 9x9 grid: C++ and Python produce similar TTT.""" + from uxsim.DTAsolvers import SolverDSO_D2D + + solver_cpp = SolverDSO_D2D(lambda: _create_dta_grid_world(seed=42, cpp=True), cpp=True) + solver_cpp.solve(max_iter=1, print_progress=False) + ttt_cpp = solver_cpp.W_sol.analyzer.total_travel_time + trips_cpp = solver_cpp.W_sol.analyzer.trip_completed + + solver_py = SolverDSO_D2D(lambda: _create_dta_grid_world(seed=42, cpp=False)) + solver_py.solve(max_iter=1, print_progress=False) + ttt_py = solver_py.W_sol.analyzer.total_travel_time + trips_py = solver_py.W_sol.analyzer.trip_completed + + assert trips_cpp == trips_py + assert ttt_cpp == pytest.approx(ttt_py, rel=0.2) + + +def test_dta_solver_due_grid_convergence_metrics_cpp(): + """SolverDUE on 9x9 grid: tracks convergence metrics (ttts, n_swaps) with cpp=True.""" + from uxsim.DTAsolvers import SolverDUE + + n_iter = 3 + solver = SolverDUE(lambda: _create_dta_grid_world(seed=42, cpp=True), cpp=True) + solver.solve(max_iter=n_iter, print_progress=False) + + assert len(solver.ttts) == n_iter + assert len(solver.n_swaps) == n_iter + assert all(t > 0 for t in solver.ttts) + assert solver.W_sol is not None + + +def test_dta_solver_dso_d2d_grid_convergence_metrics_cpp(): + """SolverDSO_D2D on 9x9 grid: tracks convergence metrics (ttts, n_swaps) with cpp=True.""" + from uxsim.DTAsolvers import SolverDSO_D2D + + n_iter = 3 + solver = SolverDSO_D2D(lambda: _create_dta_grid_world(seed=42, cpp=True), cpp=True) + solver.solve(max_iter=n_iter, print_progress=False) + + assert len(solver.ttts) == n_iter + assert len(solver.n_swaps) == n_iter + assert all(t > 0 for t in solver.ttts) + assert solver.W_sol is not None + + +def test_dta_solver_due_grid_duo_baseline_cpp(): + """DUO on 9x9 grid with cpp=True produces same TTT as Python DUO.""" + W_cpp = _create_dta_grid_world(seed=42, cpp=True) + W_cpp.exec_simulation() + W_cpp.analyzer.print_simple_stats(force_print=False) + ttt_cpp = W_cpp.analyzer.total_travel_time + + W_py = _create_dta_grid_world(seed=42, cpp=False) + W_py.exec_simulation() + W_py.analyzer.print_simple_stats(force_print=False) + ttt_py = W_py.analyzer.total_travel_time + + assert ttt_cpp > 0 + assert ttt_cpp == pytest.approx(ttt_py, rel=0.1) + + +def test_dta_solver_due_grid_benchmark_cpp(): + """Benchmark: SolverDUE on 9x9 grid with cpp=True, prints wall-clock time.""" + import time + from uxsim.DTAsolvers import SolverDUE + + t0 = time.perf_counter() + solver = SolverDUE(lambda: _create_dta_grid_world(seed=42, cpp=True), cpp=True) + solver.solve(max_iter=1, print_progress=False) + elapsed_cpp = time.perf_counter() - t0 + + t0 = time.perf_counter() + solver_py = SolverDUE(lambda: _create_dta_grid_world(seed=42, cpp=False)) + solver_py.solve(max_iter=1, print_progress=False) + elapsed_py = time.perf_counter() - t0 + + speedup = elapsed_py / elapsed_cpp if elapsed_cpp > 0 else float("inf") + print(f"\n[DTA Benchmark] SolverDUE 9x9 grid (max_iter=1):") + print(f" C++: {elapsed_cpp:.3f}s") + print(f" Python: {elapsed_py:.3f}s") + print(f" Speedup: {speedup:.2f}x") + + assert solver.W_sol.analyzer.total_travel_time > 0 + + +def test_dta_solver_dso_d2d_grid_benchmark_cpp(): + """Benchmark: SolverDSO_D2D on 9x9 grid with cpp=True, prints wall-clock time.""" + import time + from uxsim.DTAsolvers import SolverDSO_D2D + + t0 = time.perf_counter() + solver = SolverDSO_D2D(lambda: _create_dta_grid_world(seed=42, cpp=True), cpp=True) + solver.solve(max_iter=1, print_progress=False) + elapsed_cpp = time.perf_counter() - t0 + + t0 = time.perf_counter() + solver_py = SolverDSO_D2D(lambda: _create_dta_grid_world(seed=42, cpp=False)) + solver_py.solve(max_iter=1, print_progress=False) + elapsed_py = time.perf_counter() - t0 + + speedup = elapsed_py / elapsed_cpp if elapsed_cpp > 0 else float("inf") + print(f"\n[DTA Benchmark] SolverDSO_D2D 9x9 grid (max_iter=1):") + print(f" C++: {elapsed_cpp:.3f}s") + print(f" Python: {elapsed_py:.3f}s") + print(f" Speedup: {speedup:.2f}x") + + assert solver.W_sol.analyzer.total_travel_time > 0 diff --git a/uxsim/DTAsolvers/DTAsolvers.py b/uxsim/DTAsolvers/DTAsolvers.py index 99e8bff..b09f838 100644 --- a/uxsim/DTAsolvers/DTAsolvers.py +++ b/uxsim/DTAsolvers/DTAsolvers.py @@ -16,8 +16,402 @@ import warnings + +def _ensure_cpp_world(W, solver_name="CppSolver"): + """Verify that W is a CppWorld instance (has _cpp_world attribute). + + Raises ValueError if func_World did not return a World(cpp=True). + """ + if not hasattr(W, '_cpp_world'): + raise ValueError( + f"{solver_name} requires func_World to return World(cpp=True, ...). " + f"Got {type(W).__name__} without C++ engine. " + f"Add cpp=True to your World() constructor inside func_World." + ) + + +def _lightweight_finalize(W, cached_analyzer=None): + """Lightweight finalize_scenario for DTA intermediate iterations. + + Performs essential initialization (TMAX, adj_matrix, TSIZE) but reuses + a cached Analyzer object to skip expensive font loading. + + Parameters + ---------- + W : CppWorld + cached_analyzer : Analyzer or None + Analyzer from a previous iteration to reuse. If None, does full init. + + Returns + ------- + Analyzer : the (possibly reused) Analyzer object now attached to W. + """ + if W.TMAX is None: + W.TMAX = (getattr(W, '_max_demand_t', 0) // 1800 + 2) * 1800 + W._ensure_cpp_world() + W._cpp_world.initialize_adj_matrix() + + # Minimal _setup_analyzer equivalent + W.T = 0 + W.TIME = 0 + W.TSIZE = int(W.TMAX / W.DELTAT) + W.Q_AREA = defaultdict(lambda: np.zeros(int(W.TMAX / W.EULAR_DT))) + W.K_AREA = defaultdict(lambda: np.zeros(int(W.TMAX / W.EULAR_DT))) + for l in W.LINKS: + l.init_after_tmax_fix() + W.ADJ_MAT = np.zeros([len(W.NODES), len(W.NODES)]) + W.ADJ_MAT_LINKS = dict() + W.NODE_PAIR_LINKS = dict() + for link in W.LINKS: + ii = link.start_node.id + jj = link.end_node.id + W.ADJ_MAT[ii, jj] = 1 + W.ADJ_MAT_LINKS[ii, jj] = link + W.NODE_PAIR_LINKS[link.start_node.name, link.end_node.name] = link + + from ..uxsim_cpp_wrapper import CppRouteChoice + W.ROUTECHOICE = CppRouteChoice(W) + + if cached_analyzer is not None: + # Reuse existing Analyzer, just rebind to new World and reset stats + analyzer = cached_analyzer + analyzer.W = W + analyzer.average_speed = 0 + analyzer.average_speed_count = 0 + analyzer.trip_completed = 0 + analyzer.trip_all = 0 + analyzer.total_travel_time = 0 + analyzer.average_travel_time = 0 + analyzer.flag_edie_state_computed = 0 + analyzer.flag_trajectory_computed = 0 + analyzer.flag_pandas_convert = 0 + analyzer.flag_od_analysis = 0 + W.analyzer = analyzer + else: + from ..analyzer import Analyzer + W.analyzer = Analyzer(W) + + # Pre-compute congestion pricing tolls + for link in W.LINKS: + if link.congestion_pricing is not None: + tolls = [float(link.congestion_pricing(t * W.DELTAT)) for t in range(W.TSIZE)] + link._cpp_link.set_toll_timeseries(tolls) + + W.finalized = 1 + W.sim_start_time = time.time() + + return W.analyzer + + +def _build_name_id_maps(W): + """Build name<->ID mapping dicts for a CppWorld. Called once per iteration.""" + link_name_to_id = {link.name: link.id for link in W.LINKS} + link_id_to_name = [link.name for link in W.LINKS] + node_name_to_id = {node.name: node.id for node in W.NODES} + return link_name_to_id, link_id_to_name, node_name_to_id + + +def _convert_od_routes_to_ids(dict_od_to_routes, node_name_to_id, link_name_to_id): + """Convert dict_od_to_routes from str names to int IDs for C++. + + Parameters + ---------- + dict_od_to_routes : dict + {(o_name, d_name): [[link_name, ...], ...]} + node_name_to_id : dict + {node_name: node_id} + link_name_to_id : dict + {link_name: link_id} + + Returns + ------- + dict : {(node_id_o, node_id_d): [[link_id, ...], ...]} + """ + od_route_sets = {} + for (o_name, d_name), routes in dict_od_to_routes.items(): + o_id = node_name_to_id[o_name] + d_id = node_name_to_id[d_name] + od_route_sets[(o_id, d_id)] = [ + [link_name_to_id[ln] for ln in route] for route in routes + ] + return od_route_sets + + +def _convert_cpp_result_to_names(result, W, link_id_to_name): + """Convert C++ route_swap result from int IDs to str names. + + Parameters + ---------- + result : dict + C++ result with routes_specified/route_actual as list[list[int]], + cost_actual as numpy array + W : CppWorld + World object (for vehicle name list) + link_id_to_name : list + link_id -> link_name mapping + + Returns + ------- + dict with routes_specified/route_actual as dict[str, list[str]], + cost_actual as dict[str, float], and scalar metrics + """ + veh_names = list(W.VEHICLES.keys()) + rs = result['routes_specified'] + ra = result['route_actual'] + ca = result['cost_actual'] + + routes_specified_data = {} + for idx, route_ids in enumerate(rs): + if route_ids: + routes_specified_data[veh_names[idx]] = [link_id_to_name[lid] for lid in route_ids] + + route_actual = {} + for idx, route_ids in enumerate(ra): + route_actual[veh_names[idx]] = [link_id_to_name[lid] for lid in route_ids] + + cost_actual = {} + for idx, veh_name in enumerate(veh_names): + cost_actual[veh_name] = float(ca[idx]) + + return { + 'routes_specified': routes_specified_data, + 'route_actual': route_actual, + 'cost_actual': cost_actual, + 'n_swap': float(result['n_swap']), + 'potential_n_swap': float(result['potential_n_swap']), + 'total_t_gap': float(result['total_t_gap']), + } + + +def _build_enforce_routes_input(W, routes_specified_data, link_name_to_id): + """Convert routes_specified_data to list[list[int]] for C++ batch_enforce_routes. + + Parameters + ---------- + W : CppWorld + routes_specified_data : dict[str, list[str]] + {veh_name: [link_name, ...]} + link_name_to_id : dict + {link_name: link_id} + + Returns + ------- + list[list[int]] : one entry per vehicle in VEHICLES order. Empty list = skip. + """ + routes_per_vehicle = [] + for veh_name in W.VEHICLES: + if veh_name in routes_specified_data: + routes_per_vehicle.append( + [link_name_to_id[ln] for ln in routes_specified_data[veh_name]] + ) + else: + routes_per_vehicle.append([]) + return routes_per_vehicle + + +def _route_swap_due_python(W, dict_od_to_routes, dict_od_to_vehid, swap_prob, has_fixed_route_sets): + """Python fallback for DUE route swap logic. + + This function extracts the route swap loop from SolverDUE.solve() so that + it can be replaced by a C++ batch call later. The logic is identical to + the original Python loop. + + Parameters + ---------- + W : World or CppWorld + World object after simulation execution + dict_od_to_routes : dict + {(o_name, d_name): [[link_name, ...], ...]} + dict_od_to_vehid : dict + {(o_name, d_name): [veh_key, ...]} + swap_prob : float + probability of route swap per vehicle + has_fixed_route_sets : bool + whether route_sets was user-provided + + Returns + ------- + dict with keys: routes_specified, route_actual, cost_actual, + n_swap, potential_n_swap, total_t_gap + """ + # Build Route objects for this World + route_set = defaultdict(lambda: []) + for o, d in dict_od_to_vehid.keys(): + for r in dict_od_to_routes[o, d]: + route_set[o, d].append(W.defRoute(r)) + + routes_specified_data = {} + route_actual = {} + cost_actual = {} + n_swap = 0 + total_t_gap = 0 + potential_n_swap = 0 + + for key, veh in W.VEHICLES.items(): + flag_swap = random.random() < swap_prob + o = veh.orig.name + d = veh.dest.name + r, ts = veh.traveled_route() + travel_time = ts[-1] - ts[0] + + route_actual[key] = [rr.name for rr in r] + cost_actual[key] = travel_time + + if veh.state != "end": + continue + + if has_fixed_route_sets and r not in route_set[o, d]: + routes_specified_data[key] = route_set[o, d][0].links_name + continue + + flag_route_changed = False + route_changed_links = None + t_gap = 0 + + cost_current = r.actual_travel_time(ts[0]) + sum([l.get_toll(ts[i]) for i, l in enumerate(r)]) + + potential_n_swap_updated = potential_n_swap + for alt_route in route_set[o, d]: + alt_route_tts = alt_route.actual_travel_time(ts[0], return_details=True)[1] + cost_alt = alt_route.actual_travel_time(ts[0]) + sum([l.get_toll(ts[0] + sum(alt_route_tts[:i])) for i, l in enumerate(alt_route)]) + if cost_alt < cost_current: + if flag_route_changed == False or (cost_alt < cost_current): + t_gap = cost_current - cost_alt + potential_n_swap_updated = potential_n_swap + W.DELTAN + if flag_swap: + flag_route_changed = True + route_changed_links = alt_route.links_name + cost_current = cost_alt + + potential_n_swap = potential_n_swap_updated + + total_t_gap += t_gap + routes_specified_data[key] = [rr.name for rr in r] + if flag_route_changed: + n_swap += W.DELTAN + routes_specified_data[key] = route_changed_links + + return { + 'routes_specified': routes_specified_data, + 'route_actual': route_actual, + 'cost_actual': cost_actual, + 'n_swap': n_swap, + 'potential_n_swap': potential_n_swap, + 'total_t_gap': total_t_gap, + } + + +def _route_swap_dso_python(W, dict_od_to_routes, dict_od_to_vehid, swap_prob, swap_num, has_fixed_route_sets): + """Python fallback for DSO route swap logic. + + This function extracts the route swap loop from SolverDSO_D2D.solve() so + that it can be replaced by a C++ batch call later. The logic is identical + to the original Python loop, using marginal cost (private + externality). + + Parameters + ---------- + W : World or CppWorld + World object after simulation execution + dict_od_to_routes : dict + {(o_name, d_name): [[link_name, ...], ...]} + dict_od_to_vehid : dict + {(o_name, d_name): [veh_key, ...]} + swap_prob : float + probability of route swap per vehicle + swap_num : int or None + fixed number of vehicles to swap (overrides swap_prob if set) + has_fixed_route_sets : bool + whether route_sets was user-provided + + Returns + ------- + dict with keys: routes_specified, route_actual, cost_actual, + n_swap, potential_n_swap, total_t_gap + """ + # Build Route objects for this World + route_set = defaultdict(lambda: []) + for o, d in dict_od_to_vehid.keys(): + for r in dict_od_to_routes[o, d]: + route_set[o, d].append(W.defRoute(r)) + + routes_specified_data = {} + route_actual = {} + cost_actual = {} + n_swap = 0 + total_t_gap = 0 + potential_n_swap = 0 + + keys_swap = [key for key in W.VEHICLES.keys() if random.random() < swap_prob] + if swap_num is not None: + keys_swap = random.sample(list(W.VEHICLES.keys()), swap_num) + + for key, veh in W.VEHICLES.items(): + flag_swap = key in keys_swap + o = veh.orig.name + d = veh.dest.name + r, ts = veh.traveled_route() + travel_time = ts[-1] - ts[0] + + route_actual[key] = [rr.name for rr in r] + cost_actual[key] = travel_time + + if veh.state != "end": + continue + + if has_fixed_route_sets and r not in route_set[o, d]: + routes_specified_data[key] = route_set[o, d][0].links_name + continue + + flag_route_changed = False + route_changed_links = None + t_gap = 0 + + ext = estimate_congestion_externality_route(W, r, ts[0]) + private_cost = r.actual_travel_time(ts[0]) + cost_current = private_cost + ext + + potential_n_swap_updated = potential_n_swap + + for alt_route in route_set[o, d]: + ext = estimate_congestion_externality_route(W, alt_route, ts[0]) + private_cost = alt_route.actual_travel_time(ts[0]) + cost_alt = private_cost + ext + + if cost_alt < cost_current: + if flag_route_changed == False or (cost_alt < cost_current): + t_gap = cost_current - cost_alt + potential_n_swap_updated = potential_n_swap + W.DELTAN + if flag_swap: + flag_route_changed = True + route_changed_links = alt_route.links_name + cost_current = cost_alt + + potential_n_swap = potential_n_swap_updated + + total_t_gap += t_gap + routes_specified_data[key] = [rr.name for rr in r] + if flag_route_changed: + n_swap += W.DELTAN + routes_specified_data[key] = route_changed_links + + return { + 'routes_specified': routes_specified_data, + 'route_actual': route_actual, + 'cost_actual': cost_actual, + 'n_swap': n_swap, + 'potential_n_swap': potential_n_swap, + 'total_t_gap': total_t_gap, + } + + class SolverDUE: - def __init__(s, func_World): + + def __new__(cls, func_World, cpp=False): + if cpp and cls is SolverDUE: + return CppSolverDUE(func_World) + return super().__new__(cls) + + def __init__(s, func_World, cpp=False): """ Solve quasi Dynamic User Equilibrium (DUE) problem using day-to-day dynamics. @@ -25,7 +419,9 @@ def __init__(s, func_World): ---------- func_World : function function that returns a World object with nodes, links, and demand specifications - + cpp : bool + if True, return a CppSolverDUE instance that uses C++ batch operations + Notes ----- This function computes a near dynamic user equilibrium state as a steady state of day-to-day dynamical routing game. @@ -44,6 +440,8 @@ def __init__(s, func_World): - Iryo, T., Urata, J., & Kawase, R. (2024). Traffic Flow Simulator and Travel Demand Simulators for Assessing Congestion on Roads After a Major Earthquake. In APPLICATION OF HIGH-PERFORMANCE COMPUTING TO EARTHQUAKE-RELATED PROBLEMS (pp. 413-447). https://doi.org/10.1142/9781800614635_0007 - Iryo, T., Watling, D., & Hazelton, M. (2024). Estimating Markov Chain Mixing Times: Convergence Rate Towards Equilibrium of a Stochastic Process Traffic Assignment Model. Transportation Science. https://doi.org/10.1287/trsc.2024.0523 """ + if isinstance(s, CppSolverDUE): + return # CppSolverDUE.__init__ handles its own init s.func_World = func_World s.W_sol = None #final solution s.W_intermid_solution = None #latest solution in the iterative process. Can be used when a user terminates the solution algorithm @@ -89,10 +487,10 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, route_sets=None, prin ------- W : World World object with quasi DUE solution (if properly converged) - + Notes ----- - `self.W_sol` is the final solution. + `self.W_sol` is the final solution. `self.W_intermid_solution` is the latest solution in the iterative process. Can be used when a user terminates the solution algorithm. """ s.start_time = time.time() @@ -151,11 +549,11 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, route_sets=None, prin for key in W.VEHICLES: if key in routes_specified: W.VEHICLES[key].enforce_route(routes_specified[key]) - + route_set = defaultdict(lambda: []) #routes[o,d] = [Route, Route, ...] for o,d in dict_od_to_vehid.keys(): for r in dict_od_to_routes[o,d]: - route_set[o,d].append(W.defRoute(r)) + route_set[o,d].append(W.defRoute(r)) # simulation W.exec_simulation() @@ -171,7 +569,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, route_sets=None, prin # attach route choice set to W object for later re-use at different solvers like DSO-GA W.dict_od_to_routes = dict_od_to_routes - + s.W_intermid_solution = W s.dfs_link.append(W.analyzer.link_to_pandas()) @@ -189,7 +587,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, route_sets=None, prin d = veh.dest.name r, ts = veh.traveled_route() travel_time = ts[-1]-ts[0] - + route_actual[key] = [rr.name for rr in r] cost_actual[key] = travel_time @@ -206,7 +604,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, route_sets=None, prin t_gap = 0 cost_current = r.actual_travel_time(ts[0]) + sum([l.get_toll(ts[i]) for i,l in enumerate(r)]) - + potential_n_swap_updated = potential_n_swap for alt_route in route_set[o,d]: alt_route_tts = alt_route.actual_travel_time(ts[0], return_details=True)[1] @@ -219,9 +617,9 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, route_sets=None, prin flag_route_changed = True route_changed = alt_route cost_current = cost_alt - + potential_n_swap = potential_n_swap_updated - + total_t_gap += t_gap routes_specified[key] = r if flag_route_changed: @@ -239,7 +637,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, route_sets=None, prin s.n_swaps.append(n_swap) s.potential_swaps.append(potential_n_swap) s.t_gaps.append(t_gap_per_vehicle) - + s.end_time = time.time() print("DUE summary:") @@ -363,35 +761,44 @@ def plot_vehicle_stats(s, orig=None, dest=None): class SolverDSO_D2D: - def __init__(s, func_World): + + def __new__(cls, func_World, cpp=False): + if cpp and cls is SolverDSO_D2D: + return CppSolverDSO_D2D(func_World) + return super().__new__(cls) + + def __init__(s, func_World, cpp=False): """ - Solve quasi Dynamic System Optimum (DSO) problem using day-to-day dynamics. + Solve quasi Dynamic System Optimum (DSO) problem using day-to-day dynamics. Parameters ---------- func_World : function function that returns a World object with nodes, links, and demand specifications - + cpp : bool + if True, return a CppSolverDSO_D2D instance that uses C++ batch operations + Notes ----- This function computes a near DSO state as a steady state of day-to-day dynamical routing game. - + This is a modification of Iryo's DUE algorithm, in which each traveler minimize marginal travel time (i.e., private cost + externality) instead of private cost as in DUE. Since the externality (i.e., the total system cost without the ego vehicle) is costly to compute directly, it is estimated based on the queueing principle of the traffic flow model. Apart from these point, this algorithm is almost the same to `SolverDUE` and thus expected to have the desirable properties discussed in the papers cited above. - + This algorithm is a loose generalization of "DSO game" of the following paper. - Satsukawa, K., Wada, K., & Watling, D. (2022). Dynamic system optimal traffic assignment with atomic users: Convergence and stability. Transportation Research Part B: Methodological, 155, 188-209. In this paper, it is theoretically guaranteed that their DSO game algorithm converges to the global optimal DSO state. However, it may take a long time. This `SolverDSO_D2D` speeds up the solution process significantly, but it weakens the convergence guarantee. """ - + if isinstance(s, CppSolverDSO_D2D): + return s.func_World = func_World s.W_sol = None #final solution s.W_intermid_solution = None #latest solution in the iterative process. Can be used when a user terminates the solution algorithm s.dfs_link = [] - + def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_sets=None, print_progress=True): """ Solve quasi DSO problem using day-to-day dynamics. @@ -432,10 +839,10 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_ ------- W : World World object with quasi DUE solution (if properly converged) - + Notes ----- - `self.W_sol` is the final solution. + `self.W_sol` is the final solution. `self.W_intermid_solution` is the latest solution in the iterative process. Can be used when a user terminates the solution algorithm. """ s.start_time = time.time() @@ -470,7 +877,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_ dict_od_to_routes[o,d] = [] for route in routes: dict_od_to_routes[o,d].append([W_orig.get_link(l).name for l in route]) - + if print_progress: print(f"number of OD pairs: {len(dict_od_to_routes.keys())}, number of routes: {sum([len(val) for val in dict_od_to_routes.values()])}") @@ -494,7 +901,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_ for key in W.VEHICLES: if key in routes_specified: W.VEHICLES[key].enforce_route(routes_specified[key]) - + route_set = defaultdict(lambda: []) #routes[o,d] = [Route, Route, ...] for o,d in dict_od_to_vehid.keys(): for r in dict_od_to_routes[o,d]: @@ -514,7 +921,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_ # attach route choice set to W object for later re-use at different solvers like DSO-GA W.dict_od_to_routes = dict_od_to_routes - + if unfinished_trips == 0: if s.W_intermid_solution == None: s.W_intermid_solution = W @@ -541,7 +948,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_ d = veh.dest.name r, ts = veh.traveled_route() travel_time = ts[-1]-ts[0] - + route_actual[key] = [rr.name for rr in r] cost_actual[key] = travel_time @@ -562,7 +969,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_ cost_current = private_cost+ext potential_n_swap_updated = potential_n_swap - + for alt_route in route_set[o,d]: ext = estimate_congestion_externality_route(W, alt_route, ts[0]) private_cost = alt_route.actual_travel_time(ts[0]) @@ -575,8 +982,8 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_ if flag_swap: flag_route_changed = True route_changed = alt_route - cost_current = cost_alt - + cost_current = cost_alt + potential_n_swap = potential_n_swap_updated total_t_gap += t_gap @@ -596,7 +1003,7 @@ def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_ s.n_swaps.append(n_swap) s.potential_swaps.append(potential_n_swap) s.t_gaps.append(t_gap_per_vehicle) - + s.end_time = time.time() print("DSO summary:") @@ -1434,3 +1841,244 @@ def plot_vehicle_stats(s, orig=None, dest=None): plt.ylabel("travel time") plt.legend() plt.show() + + +# ============================================================================= +# C++ accelerated solvers (dispatched via __new__ in SolverDUE / SolverDSO_D2D) +# ============================================================================= + +class CppSolverDUE(SolverDUE): + """C++ accelerated DUE solver. Created via ``SolverDUE(func_World, cpp=True)``.""" + + def __init__(s, func_World, **kwargs): + s.func_World = func_World + s.W_sol = None + s.W_intermid_solution = None + s.dfs_link = [] + + def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, route_sets=None, print_progress=True): + s.start_time = time.time() + + W_orig = s.func_World() + _ensure_cpp_world(W_orig, "CppSolverDUE") + if print_progress: + W_orig.print_scenario_stats() + + dict_od_to_vehid = defaultdict(lambda: []) + for key, veh in W_orig.VEHICLES.items(): + dict_od_to_vehid[veh.orig.name, veh.dest.name].append(key) + + if W_orig.finalized == False: + W_orig.finalize_scenario() + dict_od_to_routes = enumerate_k_random_routes(W_orig, k=n_routes_per_od) + + if route_sets is not None: + dict_od_to_routes = {} + for key, routes in route_sets.items(): + o = W_orig.get_node(key[0]).name + d = W_orig.get_node(key[1]).name + dict_od_to_routes[o, d] = [[W_orig.get_link(l).name for l in route] for route in routes] + + if print_progress: + print(f"number of OD pairs: {len(dict_od_to_routes.keys())}, number of routes: {sum([len(val) for val in dict_od_to_routes.values()])}") + + s.ttts = [] + s.n_swaps = [] + s.potential_swaps = [] + s.t_gaps = [] + s.route_log = [] + s.cost_log = [] + + od_route_sets_ids = None + link_name_to_id = None + link_id_to_name = None + cached_analyzer = None + + print("solving DUE...") + for i in range(max_iter): + W = s.func_World() + is_last_iter = (i == max_iter - 1) + + if od_route_sets_ids is None: + link_name_to_id, link_id_to_name, node_name_to_id = _build_name_id_maps(W) + od_route_sets_ids = _convert_od_routes_to_ids(dict_od_to_routes, node_name_to_id, link_name_to_id) + + cached_analyzer = _lightweight_finalize(W, cached_analyzer) + + # Skip expensive log cache building on intermediate iterations; + # CppWorld log mode is set at __init__ time, so vehicle_logging_timestep_interval + # cannot be changed after construction — _skip_log_on_terminate handles this instead. + if not is_last_iter: + W._skip_log_on_terminate = True + + if i != 0: + enforce_input = _build_enforce_routes_input(W, routes_specified_data, link_name_to_id) + W._cpp_world.batch_enforce_routes(enforce_input) + + W.exec_simulation() + W.analyzer.print_simple_stats() + + unfinished_trips = W.analyzer.trip_all - W.analyzer.trip_completed + if unfinished_trips > 0: + warnings.warn(f"Warning: {unfinished_trips} / {W.analyzer.trip_all} vehicles have not finished their trips. The DUE solver assumes that all vehicles finish their trips during the simulation duration. Consider increasing the simulation time limit or checking the network configuration.", UserWarning) + + W.dict_od_to_routes = dict_od_to_routes + s.W_intermid_solution = W + if is_last_iter: + s.dfs_link.append(W.analyzer.link_to_pandas()) + + rng_seed = random.randint(0, 2**31 - 1) + cpp_result = W._cpp_world.route_swap_due( + od_route_sets_ids, swap_prob, route_sets is not None, rng_seed + ) + result = _convert_cpp_result_to_names(cpp_result, W, link_id_to_name) + routes_specified_data = result['routes_specified'] + route_actual = result['route_actual'] + cost_actual = result['cost_actual'] + n_swap = result['n_swap'] + potential_n_swap = result['potential_n_swap'] + total_t_gap = result['total_t_gap'] + + t_gap_per_vehicle = total_t_gap / len(W.VEHICLES) + if print_progress: + print(f' iter {i}: time gap: {t_gap_per_vehicle:.1f}, potential route change: {potential_n_swap}, route change: {n_swap}, total travel time: {W.analyzer.total_travel_time: .1f}, delay ratio: {W.analyzer.average_delay/W.analyzer.average_travel_time: .3f}') + + s.route_log.append(route_actual) + s.cost_log.append(cost_actual) + s.ttts.append(int(W.analyzer.total_travel_time)) + s.n_swaps.append(n_swap) + s.potential_swaps.append(potential_n_swap) + s.t_gaps.append(t_gap_per_vehicle) + + s.end_time = time.time() + + print("DUE summary:") + last_iters = int(max_iter / 4) + print(f" total travel time: initial {s.ttts[0]:.1f} -> average of last {last_iters} iters {np.average(s.ttts[-last_iters:]):.1f}") + print(f" number of potential route changes: initial {s.potential_swaps[0]:.1f} -> average of last {last_iters} iters {np.average(s.potential_swaps[-last_iters:]):.1f}") + print(f" route travel cost gap: initial {s.t_gaps[0]:.1f} -> average of last {last_iters} iters {np.average(s.t_gaps[-last_iters:]):.1f}") + print(f" computation time: {s.end_time - s.start_time:.1f} seconds") + + s.W_sol = W + return s.W_sol + + +class CppSolverDSO_D2D(SolverDSO_D2D): + """C++ accelerated DSO solver. Created via ``SolverDSO_D2D(func_World, cpp=True)``.""" + + def __init__(s, func_World, **kwargs): + s.func_World = func_World + s.W_sol = None + s.W_intermid_solution = None + s.dfs_link = [] + + def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, swap_num=None, route_sets=None, print_progress=True): + s.start_time = time.time() + + W_orig = s.func_World() + _ensure_cpp_world(W_orig, "CppSolverDSO_D2D") + if print_progress: + W_orig.print_scenario_stats() + + dict_od_to_vehid = defaultdict(lambda: []) + for key, veh in W_orig.VEHICLES.items(): + dict_od_to_vehid[veh.orig.name, veh.dest.name].append(key) + + if W_orig.finalized == False: + W_orig.finalize_scenario() + dict_od_to_routes = enumerate_k_random_routes(W_orig, k=n_routes_per_od) + + if route_sets is not None: + dict_od_to_routes = {} + for key, routes in route_sets.items(): + o = W_orig.get_node(key[0]).name + d = W_orig.get_node(key[1]).name + dict_od_to_routes[o, d] = [[W_orig.get_link(l).name for l in route] for route in routes] + + if print_progress: + print(f"number of OD pairs: {len(dict_od_to_routes.keys())}, number of routes: {sum([len(val) for val in dict_od_to_routes.values()])}") + + s.ttts = [] + s.n_swaps = [] + s.potential_swaps = [] + s.t_gaps = [] + s.route_log = [] + s.cost_log = [] + + od_route_sets_ids = None + link_name_to_id = None + link_id_to_name = None + cached_analyzer = None + + print("solving DSO...") + for i in range(max_iter): + W = s.func_World() + is_last_iter = (i == max_iter - 1) + + if od_route_sets_ids is None: + link_name_to_id, link_id_to_name, node_name_to_id = _build_name_id_maps(W) + od_route_sets_ids = _convert_od_routes_to_ids(dict_od_to_routes, node_name_to_id, link_name_to_id) + + cached_analyzer = _lightweight_finalize(W, cached_analyzer) + + # CppWorld log mode is fixed at __init__; use _skip_log_on_terminate instead + if not is_last_iter: + W._skip_log_on_terminate = True + + if i != 0: + enforce_input = _build_enforce_routes_input(W, routes_specified_data, link_name_to_id) + W._cpp_world.batch_enforce_routes(enforce_input) + + W.exec_simulation() + W.analyzer.print_simple_stats() + + unfinished_trips = W.analyzer.trip_all - W.analyzer.trip_completed + if unfinished_trips > 0: + warnings.warn(f"Warning: {unfinished_trips} / {W.analyzer.trip_all} vehicles have not finished their trips. The DSO solver assumes that all vehicles finish their trips during the simulation duration. Consider increasing the simulation time limit or checking the network configuration.", UserWarning) + + W.dict_od_to_routes = dict_od_to_routes + + if unfinished_trips == 0: + if s.W_intermid_solution is None: + s.W_intermid_solution = W + elif W.analyzer.average_travel_time < s.W_intermid_solution.analyzer.average_travel_time: + s.W_intermid_solution = W + + if is_last_iter: + s.dfs_link.append(W.analyzer.link_to_pandas()) + + rng_seed = random.randint(0, 2**31 - 1) + cpp_swap_num = swap_num if swap_num is not None else -1 + cpp_result = W._cpp_world.route_swap_dso( + od_route_sets_ids, swap_prob, cpp_swap_num, route_sets is not None, rng_seed + ) + result = _convert_cpp_result_to_names(cpp_result, W, link_id_to_name) + routes_specified_data = result['routes_specified'] + route_actual = result['route_actual'] + cost_actual = result['cost_actual'] + n_swap = result['n_swap'] + potential_n_swap = result['potential_n_swap'] + total_t_gap = result['total_t_gap'] + + t_gap_per_vehicle = total_t_gap / len(W.VEHICLES) + if print_progress: + print(f' iter {i}: marginal time gap: {t_gap_per_vehicle:.1f}, potential route change: {potential_n_swap}, route change: {n_swap}, total travel time: {W.analyzer.total_travel_time: .1f}, delay ratio: {W.analyzer.average_delay/W.analyzer.average_travel_time: .3f}') + + s.route_log.append(route_actual) + s.cost_log.append(cost_actual) + s.ttts.append(int(W.analyzer.total_travel_time)) + s.n_swaps.append(n_swap) + s.potential_swaps.append(potential_n_swap) + s.t_gaps.append(t_gap_per_vehicle) + + s.end_time = time.time() + + print("DSO summary:") + last_iters = int(max_iter / 4) + print(f" total travel time: initial {s.ttts[0]:.1f} -> minimum {np.min(s.ttts):.1f}") + print(f" number of potential route changes: initial {s.potential_swaps[0]:.1f} -> minimum {np.min(s.potential_swaps):.1f}") + print(f" route marginal travel time gap: initial {s.t_gaps[0]:.1f} -> minimum {np.min(s.t_gaps):.1f}") + print(f" computation time: {s.end_time - s.start_time:.1f} seconds") + + s.W_sol = s.W_intermid_solution + return s.W_sol diff --git a/uxsim/trafficpp/bindings.cpp b/uxsim/trafficpp/bindings.cpp index c455e5c..e2d83fe 100644 --- a/uxsim/trafficpp/bindings.cpp +++ b/uxsim/trafficpp/bindings.cpp @@ -16,6 +16,7 @@ #include #include "traffi.cpp" +#include "dta_solver.cpp" namespace nb = nanobind; @@ -523,6 +524,112 @@ NB_MODULE(uxsim_cpp, m) { return d; }, "Build enter_log data for all vehicles. Returns dict with link_id/time/vehicle_index arrays.") + // ----- DTA Solver batch operations (from dta_solver.h) ----- + .def("batch_enforce_routes", [](World &w, nb::list py_routes) { + size_t nv = nb::len(py_routes); + std::vector> routes(nv); + for (size_t i = 0; i < nv; i++) { + nb::list r = nb::cast(py_routes[i]); + size_t nr = nb::len(r); + routes[i].reserve(nr); + for (size_t j = 0; j < nr; j++) { + routes[i].push_back(nb::cast(r[j])); + } + } + dta_batch_enforce_routes(&w, routes); + }, + nb::arg("routes_per_vehicle"), + "Batch enforce routes for all vehicles. routes_per_vehicle[i] = list of link IDs (empty = skip)") + .def("route_swap_due", [](World &w, nb::dict py_od_routes, double swap_prob, + bool has_external_route_sets, unsigned int rng_seed) -> nb::dict { + std::map, std::vector>> od_route_sets; + for (auto [key, val] : py_od_routes) { + nb::tuple k = nb::cast(key); + int o = nb::cast(k[0]); + int d = nb::cast(k[1]); + nb::list routes_list = nb::cast(val); + std::vector> routes; + for (size_t ri = 0; ri < nb::len(routes_list); ri++) { + nb::list route = nb::cast(routes_list[ri]); + std::vector r; + for (size_t li = 0; li < nb::len(route); li++) { + r.push_back(nb::cast(route[li])); + } + routes.push_back(std::move(r)); + } + od_route_sets[{o, d}] = std::move(routes); + } + auto result = dta_route_swap_due(&w, od_route_sets, swap_prob, has_external_route_sets, rng_seed); + nb::dict out; + nb::list rs_list; + for (auto &r : result.routes_specified) { + nb::list rl; + for (int lid : r) rl.append(lid); + rs_list.append(rl); + } + out["routes_specified"] = rs_list; + nb::list ra_list; + for (auto &r : result.route_actual) { + nb::list rl; + for (int lid : r) rl.append(lid); + ra_list.append(rl); + } + out["route_actual"] = ra_list; + out["cost_actual"] = make_numpy_move(std::move(result.cost_actual)); + out["n_swap"] = result.n_swap; + out["potential_n_swap"] = result.potential_n_swap; + out["total_t_gap"] = result.total_t_gap; + return out; + }, + nb::arg("od_route_sets"), nb::arg("swap_prob"), + nb::arg("has_external_route_sets"), nb::arg("rng_seed"), + "DUE route swap: compute next-iteration routes for all vehicles in batch") + .def("route_swap_dso", [](World &w, nb::dict py_od_routes, double swap_prob, + int swap_num, bool has_external_route_sets, + unsigned int rng_seed) -> nb::dict { + std::map, std::vector>> od_route_sets; + for (auto [key, val] : py_od_routes) { + nb::tuple k = nb::cast(key); + int o = nb::cast(k[0]); + int d = nb::cast(k[1]); + nb::list routes_list = nb::cast(val); + std::vector> routes; + for (size_t ri = 0; ri < nb::len(routes_list); ri++) { + nb::list route = nb::cast(routes_list[ri]); + std::vector r; + for (size_t li = 0; li < nb::len(route); li++) { + r.push_back(nb::cast(route[li])); + } + routes.push_back(std::move(r)); + } + od_route_sets[{o, d}] = std::move(routes); + } + auto result = dta_route_swap_dso(&w, od_route_sets, swap_prob, swap_num, + has_external_route_sets, rng_seed); + nb::dict out; + nb::list rs_list; + for (auto &r : result.routes_specified) { + nb::list rl; + for (int lid : r) rl.append(lid); + rs_list.append(rl); + } + out["routes_specified"] = rs_list; + nb::list ra_list; + for (auto &r : result.route_actual) { + nb::list rl; + for (int lid : r) rl.append(lid); + ra_list.append(rl); + } + out["route_actual"] = ra_list; + out["cost_actual"] = make_numpy_move(std::move(result.cost_actual)); + out["n_swap"] = result.n_swap; + out["potential_n_swap"] = result.potential_n_swap; + out["total_t_gap"] = result.total_t_gap; + return out; + }, + nb::arg("od_route_sets"), nb::arg("swap_prob"), nb::arg("swap_num"), + nb::arg("has_external_route_sets"), nb::arg("rng_seed"), + "DSO route swap: compute next-iteration routes using marginal cost (private + externality)") ; // diff --git a/uxsim/trafficpp/dta_solver.cpp b/uxsim/trafficpp/dta_solver.cpp new file mode 100644 index 0000000..47d97c0 --- /dev/null +++ b/uxsim/trafficpp/dta_solver.cpp @@ -0,0 +1,346 @@ +// clang-format off + +// DTA solver batch operations — implementation. +// Separate from traffi.cpp to keep normal simulation code untouched. + +#include "dta_solver.h" + +// ----------------------------------------------------------------------- +// Traveled route extraction +// ----------------------------------------------------------------------- + +TraveledRouteInfo dta_get_traveled_route(const Vehicle *veh) { + TraveledRouteInfo info; + info.arrival_time = -1.0; + + int prev_lid = -999; + for (size_t i = 0; i < veh->log_size; i++){ + int lid = veh->log_link[i]; + if (lid >= 0 && lid != prev_lid){ + info.link_ids.push_back(lid); + info.entry_times.push_back(veh->log_t[i]); + prev_lid = lid; + } + } + + if (veh->state == vsEND && veh->arrival_time >= 0){ + info.arrival_time = veh->arrival_time; + } + return info; +} + +// ----------------------------------------------------------------------- +// Route cost computation +// ----------------------------------------------------------------------- + +double dta_compute_route_travel_time(const World *w, const std::vector &route_link_ids, + double departure_time, std::vector *details) { + double tt_total = 0.0; + double t = departure_time; + for (int lid : route_link_ids){ + Link *ln = w->links[lid]; + double ltt = dta_get_actual_travel_time(ln, t); + tt_total += ltt; + t += ltt; + if (details) details->push_back(ltt); + } + return tt_total; +} + +double dta_compute_route_cost_due(const World *w, const std::vector &route_link_ids, + double departure_time) { + double tt_total = 0.0; + double toll_total = 0.0; + double t = departure_time; + for (int lid : route_link_ids){ + Link *ln = w->links[lid]; + double ltt = dta_get_actual_travel_time(ln, t); + toll_total += dta_get_toll(ln, t); + tt_total += ltt; + t += ltt; + } + return tt_total + toll_total; +} + +double dta_compute_route_cost_dso(World *w, const std::vector &route_link_ids, + double departure_time) { + double tt_total = 0.0; + double ext_total = 0.0; + double t = departure_time; + for (int lid : route_link_ids){ + Link *ln = w->links[lid]; + double ltt = dta_get_actual_travel_time(ln, t); + int ts = (int)(t / w->delta_t); + ext_total += ln->estimate_congestion_externality(ts); + tt_total += ltt; + t += ltt; + } + return tt_total + ext_total; +} + +// ----------------------------------------------------------------------- +// Batch enforce routes +// ----------------------------------------------------------------------- + +void dta_batch_enforce_routes(World *w, const std::vector> &routes_per_vehicle) { + size_t n = std::min(routes_per_vehicle.size(), w->vehicles.size()); + for (size_t i = 0; i < n; i++){ + const auto &route_ids = routes_per_vehicle[i]; + if (route_ids.empty()) continue; + // Convert link IDs to Link pointers + std::vector route; + route.reserve(route_ids.size()); + for (int lid : route_ids){ + if (lid >= 0 && lid < (int)w->links.size()){ + route.push_back(w->links[lid]); + } + } + w->vehicles[i]->enforce_route(route); + } +} + +// ----------------------------------------------------------------------- +// DUE route swap +// ----------------------------------------------------------------------- + +RouteSwapResult dta_route_swap_due( + World *w, + const std::map, std::vector>> &od_route_sets, + double swap_prob, + bool has_external_route_sets, + unsigned int rng_seed) +{ + size_t nv = w->vehicles.size(); + RouteSwapResult result; + result.routes_specified.resize(nv); + result.route_actual.resize(nv); + result.cost_actual.resize(nv, 0.0); + result.n_swap = 0; + result.potential_n_swap = 0; + result.total_t_gap = 0; + + std::mt19937 local_rng(rng_seed); + std::uniform_real_distribution udist(0.0, 1.0); + + for (size_t vi = 0; vi < nv; vi++){ + Vehicle *veh = w->vehicles[vi]; + int orig_id = veh->orig->id; + int dest_id = veh->dest->id; + + // Get traveled route + auto traveled = dta_get_traveled_route(veh); + result.route_actual[vi] = traveled.link_ids; + + // Compute travel time (arrival - first entry) + if (!traveled.entry_times.empty() && traveled.arrival_time >= 0){ + result.cost_actual[vi] = traveled.arrival_time - traveled.entry_times[0]; + } else { + result.cost_actual[vi] = -1; + } + + // Skip vehicles that didn't finish their trip (leave empty = no enforce_route) + if (veh->state != vsEND){ + continue; + } + + // Look up route set for this OD pair + auto it = od_route_sets.find({orig_id, dest_id}); + if (it == od_route_sets.end()){ + result.routes_specified[vi] = traveled.link_ids; + continue; + } + const auto &route_set = it->second; + + // Check if external route_sets and traveled route is not in set + if (has_external_route_sets){ + bool found = false; + for (const auto &rs : route_set){ + if (rs == traveled.link_ids){ + found = true; + break; + } + } + if (!found && !route_set.empty()){ + result.routes_specified[vi] = route_set[0]; + continue; + } + } + + // Determine if this vehicle is a swap candidate + bool flag_swap = (udist(local_rng) < swap_prob); + + // Compute current route cost (DUE: travel_time + tolls at actual entry times) + double departure_t = traveled.entry_times.empty() ? veh->departure_time : traveled.entry_times[0]; + double cost_current = 0.0; + { + double t = departure_t; + for (size_t j = 0; j < traveled.link_ids.size(); j++){ + Link *ln = w->links[traveled.link_ids[j]]; + double ltt = dta_get_actual_travel_time(ln, t); + // Toll at actual entry time for current route + double toll_t = (j < traveled.entry_times.size()) ? traveled.entry_times[j] : t; + cost_current += dta_get_toll(ln, toll_t); + cost_current += ltt; + t += ltt; + } + } + + bool flag_route_changed = false; + int route_changed_idx = -1; + double t_gap = 0.0; + double potential_n_swap_delta = 0.0; + + for (size_t ri = 0; ri < route_set.size(); ri++){ + double cost_alt = dta_compute_route_cost_due(w, route_set[ri], departure_t); + if (cost_alt < cost_current){ + if (!flag_route_changed || cost_alt < cost_current){ + t_gap = cost_current - cost_alt; + potential_n_swap_delta = w->delta_n; + if (flag_swap){ + flag_route_changed = true; + route_changed_idx = (int)ri; + cost_current = cost_alt; + } + } + } + } + + result.potential_n_swap += potential_n_swap_delta; + result.total_t_gap += t_gap; + + if (flag_route_changed && route_changed_idx >= 0){ + result.routes_specified[vi] = route_set[route_changed_idx]; + result.n_swap += w->delta_n; + } else { + result.routes_specified[vi] = traveled.link_ids; + } + } + + return result; +} + +// ----------------------------------------------------------------------- +// DSO route swap +// ----------------------------------------------------------------------- + +RouteSwapResult dta_route_swap_dso( + World *w, + const std::map, std::vector>> &od_route_sets, + double swap_prob, + int swap_num, + bool has_external_route_sets, + unsigned int rng_seed) +{ + size_t nv = w->vehicles.size(); + RouteSwapResult result; + result.routes_specified.resize(nv); + result.route_actual.resize(nv); + result.cost_actual.resize(nv, 0.0); + result.n_swap = 0; + result.potential_n_swap = 0; + result.total_t_gap = 0; + + std::mt19937 local_rng(rng_seed); + std::uniform_real_distribution udist(0.0, 1.0); + + // Determine swap candidates + std::vector is_swap_candidate(nv, false); + if (swap_num >= 0){ + // Fixed number of swaps: sample swap_num indices + std::vector indices(nv); + for (size_t i = 0; i < nv; i++) indices[i] = (int)i; + std::shuffle(indices.begin(), indices.end(), local_rng); + int actual_swap_num = std::min(swap_num, (int)nv); + for (int i = 0; i < actual_swap_num; i++){ + is_swap_candidate[indices[i]] = true; + } + } else { + // Probabilistic swap + for (size_t i = 0; i < nv; i++){ + is_swap_candidate[i] = (udist(local_rng) < swap_prob); + } + } + + for (size_t vi = 0; vi < nv; vi++){ + Vehicle *veh = w->vehicles[vi]; + int orig_id = veh->orig->id; + int dest_id = veh->dest->id; + + // Get traveled route + auto traveled = dta_get_traveled_route(veh); + result.route_actual[vi] = traveled.link_ids; + + // Compute travel time + if (!traveled.entry_times.empty() && traveled.arrival_time >= 0){ + result.cost_actual[vi] = traveled.arrival_time - traveled.entry_times[0]; + } else { + result.cost_actual[vi] = -1; + } + + // Skip vehicles that didn't finish their trip (leave empty = no enforce_route) + if (veh->state != vsEND){ + continue; + } + + // Look up route set for this OD pair + auto it = od_route_sets.find({orig_id, dest_id}); + if (it == od_route_sets.end()){ + result.routes_specified[vi] = traveled.link_ids; + continue; + } + const auto &route_set = it->second; + + // Check if external route_sets and traveled route is not in set + if (has_external_route_sets){ + bool found = false; + for (const auto &rs : route_set){ + if (rs == traveled.link_ids){ + found = true; + break; + } + } + if (!found && !route_set.empty()){ + result.routes_specified[vi] = route_set[0]; + continue; + } + } + + bool flag_swap = is_swap_candidate[vi]; + + // Compute current route cost (DSO: travel_time + congestion externality) + double departure_t = traveled.entry_times.empty() ? veh->departure_time : traveled.entry_times[0]; + double cost_current = dta_compute_route_cost_dso(w, traveled.link_ids, departure_t); + + bool flag_route_changed = false; + int route_changed_idx = -1; + double t_gap = 0.0; + double potential_n_swap_delta = 0.0; + + for (size_t ri = 0; ri < route_set.size(); ri++){ + double cost_alt = dta_compute_route_cost_dso(w, route_set[ri], departure_t); + if (cost_alt < cost_current){ + if (!flag_route_changed || cost_alt < cost_current){ + t_gap = cost_current - cost_alt; + potential_n_swap_delta = w->delta_n; + if (flag_swap){ + flag_route_changed = true; + route_changed_idx = (int)ri; + cost_current = cost_alt; + } + } + } + } + + result.potential_n_swap += potential_n_swap_delta; + result.total_t_gap += t_gap; + + if (flag_route_changed && route_changed_idx >= 0){ + result.routes_specified[vi] = route_set[route_changed_idx]; + result.n_swap += w->delta_n; + } else { + result.routes_specified[vi] = traveled.link_ids; + } + } + + return result; +} diff --git a/uxsim/trafficpp/dta_solver.h b/uxsim/trafficpp/dta_solver.h new file mode 100644 index 0000000..21aa9c4 --- /dev/null +++ b/uxsim/trafficpp/dta_solver.h @@ -0,0 +1,113 @@ +// clang-format off +#pragma once + +// DTA (Dynamic Traffic Assignment) solver batch operations for C++ engine. +// Separate from traffi.h/traffi.cpp to avoid impacting normal simulation. +// Uses World/Vehicle/Link/Node via their public interfaces only. + +#include "traffi.h" + +#include +#include +#include +#include + +// ----------------------------------------------------------------------- +// Helper functions (link-level) +// ----------------------------------------------------------------------- + +// Get toll value at time t (from toll_timeseries or route_choice_penalty) +inline double dta_get_toll(const Link *ln, double t) { + if (!ln->toll_timeseries.empty()){ + int ts = (int)(t / ln->w->delta_t); + if (ts < 0) ts = 0; + if (ts >= (int)ln->toll_timeseries.size()) ts = (int)ln->toll_timeseries.size() - 1; + return ln->toll_timeseries[ts]; + } + return ln->route_choice_penalty; +} + +// Get actual travel time at time t from traveltime_real +inline double dta_get_actual_travel_time(const Link *ln, double t) { + int n = (int)ln->traveltime_real.size(); + if (n == 0) return ln->length / ln->vmax; + int idx = (int)(t / ln->w->delta_t); + if (idx < 0) idx = 0; + if (idx >= n) idx = n - 1; + return ln->traveltime_real[idx]; +} + +// ----------------------------------------------------------------------- +// Traveled route extraction +// ----------------------------------------------------------------------- + +struct TraveledRouteInfo { + std::vector link_ids; // sequence of link IDs traveled + std::vector entry_times; // time of entering each link + double arrival_time; // arrival time at destination (-1 if not arrived) +}; + +// Extract traveled route from vehicle log data +TraveledRouteInfo dta_get_traveled_route(const Vehicle *veh); + +// ----------------------------------------------------------------------- +// Route cost computation +// ----------------------------------------------------------------------- + +// Compute actual travel time along a route from departure_time. +// If details is non-null, fills per-link travel times. +double dta_compute_route_travel_time(const World *w, const std::vector &route_link_ids, + double departure_time, std::vector *details = nullptr); + +// DUE cost: actual_travel_time + sum of tolls +double dta_compute_route_cost_due(const World *w, const std::vector &route_link_ids, + double departure_time); + +// DSO cost: actual_travel_time + congestion externality +double dta_compute_route_cost_dso(World *w, const std::vector &route_link_ids, + double departure_time); + +// ----------------------------------------------------------------------- +// Batch operations +// ----------------------------------------------------------------------- + +// Batch enforce routes for all vehicles. +// routes_per_vehicle[i] = route (link_ids) for vehicle i. Empty = skip. +void dta_batch_enforce_routes(World *w, const std::vector> &routes_per_vehicle); + +// ----------------------------------------------------------------------- +// Route swap result +// ----------------------------------------------------------------------- + +struct RouteSwapResult { + // Per-vehicle: route for next iteration (link_ids). Empty = keep traveled route. + std::vector> routes_specified; + // Per-vehicle: actually traveled route (link_ids) + std::vector> route_actual; + // Per-vehicle: travel cost (travel_time = arrival - departure) + std::vector cost_actual; + // Aggregate statistics + double n_swap; + double potential_n_swap; + double total_t_gap; +}; + +// DUE route swap: compute next-iteration routes for all vehicles. +// od_route_sets: (orig_node_id, dest_node_id) -> list of routes (each = vector of link_ids) +// has_external_route_sets: if true, vehicles whose traveled route is not in set get forced to first route +RouteSwapResult dta_route_swap_due( + World *w, + const std::map, std::vector>> &od_route_sets, + double swap_prob, + bool has_external_route_sets, + unsigned int rng_seed); + +// DSO route swap: uses marginal cost (private + externality). +// swap_num: if >= 0, exactly this many vehicles are swap candidates (overrides swap_prob) +RouteSwapResult dta_route_swap_dso( + World *w, + const std::map, std::vector>> &od_route_sets, + double swap_prob, + int swap_num, + bool has_external_route_sets, + unsigned int rng_seed); diff --git a/uxsim/utils.py b/uxsim/utils.py index 2526e48..60e82f0 100644 --- a/uxsim/utils.py +++ b/uxsim/utils.py @@ -89,10 +89,12 @@ def printtry(*args, **kwargs): print("EXCEPTION:", e) +_cached_matplotlib_font = None + def get_font_for_matplotlib(fname=None): """ Get a multi-language (currently Japanese only) font for matplotlib. - TODO: check if it works on different OS. It is only checked on Japanese Windows. + TODO: check if it works on different OS. It is only checked on Japanese Windows. TODO: explore if it can be extended for other languages. Returns @@ -100,6 +102,10 @@ def get_font_for_matplotlib(fname=None): str The name of a Japanese font that is available on the system. If no Japanese font is found, "monospace" is returned. """ + global _cached_matplotlib_font + if fname is None and _cached_matplotlib_font is not None: + return _cached_matplotlib_font + from matplotlib import font_manager font_list = font_manager.findSystemFonts() @@ -112,7 +118,9 @@ def get_font_for_matplotlib(fname=None): font = "MS Gothic" else: font = "monospace" - + + if fname is None: + _cached_matplotlib_font = font return font def print_columns(*lists): diff --git a/uxsim/uxsim_cpp_wrapper.py b/uxsim/uxsim_cpp_wrapper.py index 432ccb3..2cbea90 100644 --- a/uxsim/uxsim_cpp_wrapper.py +++ b/uxsim/uxsim_cpp_wrapper.py @@ -782,6 +782,8 @@ def __init__(self, name="", deltan=5, reaction_time=1, if vehicle_logging_timestep_interval not in (0, 1, -1): warnings.warn("vehicle_logging_timestep_interval is not 0, 1, or -1. C++ mode only supports 0 (no logging) and 1 (full logging). The value will be treated as 1 (logging enabled).", stacklevel=2) self._simulation_done = False + self._skip_log_on_terminate = False + self._skip_analyzer_setup = False self._cpp_veh_registered_count = 0 def _ensure_cpp_world(self): @@ -981,6 +983,25 @@ def _setup_analyzer(self): self.T = 0 self.TIME = 0 self.TSIZE = int(self.TMAX / self.DELTAT) + + if self._skip_analyzer_setup: + # Lightweight path for DTA intermediate iterations: + # Skip edie matrices (init_after_tmax_fix) and area containers + # that are only needed for visualization/compute_edie_state. + # ADJ_MAT and Analyzer are still needed for basic_analysis/link_to_pandas. + self.ADJ_MAT = np.zeros([len(self.NODES), len(self.NODES)]) + self.ADJ_MAT_LINKS = dict() + self.NODE_PAIR_LINKS = dict() + for link in self.LINKS: + i = link.start_node.id + j = link.end_node.id + self.ADJ_MAT[i, j] = 1 + self.ADJ_MAT_LINKS[i, j] = link + self.NODE_PAIR_LINKS[link.start_node.name, link.end_node.name] = link + self.ROUTECHOICE = CppRouteChoice(self) + self.analyzer = Analyzer(self) + return + self.Q_AREA = ddict(lambda: np.zeros(int(self.TMAX / self.EULAR_DT))) self.K_AREA = ddict(lambda: np.zeros(int(self.TMAX / self.EULAR_DT))) for l in self.LINKS: @@ -1180,8 +1201,9 @@ def _build_all_vehicle_log_caches(self): def simulation_terminated(self): self.print(" simulation finished") self._simulation_done = True - self._build_vehicles_enter_log() - self._build_all_vehicle_log_caches() + if not self._skip_log_on_terminate: + self._build_vehicles_enter_log() + self._build_all_vehicle_log_caches() self.analyzer.basic_analysis() def _build_vehicles_enter_log(self):