Update routing example; Add some functions

demos_and_examples/example_22en_routing_optimization.py

@@ -0,0 +1,165 @@
+import random
+from deap import base, creator, tools, algorithms
+from pylab import *
+from uxsim import *
+from uxsim import Utilities
+
+##############################################################
+# Define UXsim World
+W = World(
+    name="", 
+    deltan=5, 
+    tmax=1200,
+    print_mode=0, save_mode=1, show_mode=1, #print is off, otherwise it will be very verbose during the genetic algorithm
+    random_seed=None,
+    duo_update_time=60
+)
+
+## generate grid network
+Utilities.generate_grid_network(W, 3, 3, length=1000)
+W.show_network()
+
+## set demand
+od_pairs = [
+    ("n(0, 0)", "n(2, 2)"),
+    ("n(2, 0)", "n(0, 2)"),
+    ("n(0, 2)", "n(2, 0)"),
+    ("n(2, 2)", "n(0, 0)"),
+]
+for od_pair in od_pairs:
+    W.adddemand(od_pair[0], od_pair[1], 0, 500, 1)
+
+W_orig = W.copy()
+
+##############################################################
+# Compute DUO as a reference
+print("Deriving DUO")
+W.exec_simulation()
+display(W.analyzer.basic_to_pandas())
+
+W_duo= W.copy()
+
+##############################################################
+# enumerate some routes between each OD pair
+routes = {}
+n_routes = 6
+for od_pair in od_pairs:
+    routes[od_pair] = Utilities.enumerate_k_shortest_routes(W, od_pair[0], od_pair[1], n_routes)
+
+for key in routes:
+    for route in routes[key]:
+        print(key, route)
+
+##############################################################
+# Prepare genetic algorithm using DEAP
+# evaluate fitness by total travel time
+def evaluate_by_total_travel_time(W):
+    W.exec_simulation()
+    print(W.analyzer.total_travel_time, end=" ")
+    return - W.analyzer.total_travel_time,
+
+# specify routing based on genome
+def specify_routes(W, genome):
+    veh_list = list(W.VEHICLES.values())
+    for i, value in enumerate(genome):
+        veh = veh_list[i]
+        veh.set_links_prefer(routes[(veh.orig.name, veh.dest.name)][value])
+
+creator.create("FitnessMax", base.Fitness, weights=(1.0,))
+creator.create("Individual", list, fitness=creator.FitnessMax)
+
+toolbox = base.Toolbox()
+
+# Initialize the individual
+n_gene = len(W.VEHICLES)
+toolbox.register("attr_gene", random.randint, 0, 5)
+toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_gene, n=n_gene)
+toolbox.register("population", tools.initRepeat, list, toolbox.individual)
+
+# Define the evaluation, crossover, and mutation functions
+toolbox.register("evaluate", evaluate_by_total_travel_time)
+toolbox.register("mate", tools.cxTwoPoint)
+toolbox.register("mutate", tools.mutFlipBit, indpb=0.1)
+toolbox.register("select", tools.selTournament, tournsize=3)
+
+
+##############################################################
+# Execute genetic algorithm
+print("Deriving DSO using genetic algorithm")
+NPOP = 30
+CXPB, MUTPB = 0.5, 0.2
+NGEN = 30
+
+# Initial population
+pop = toolbox.population(n=NPOP) 
+for ind in pop:
+    W = W_orig.copy()
+    specify_routes(W, ind)
+    ind.W = W
+fitnesses = list(map(toolbox.evaluate, [ind.W for ind in pop]))
+for ind, fit in zip(pop, fitnesses):
+    ind.fitness.values = fit
+
+for g in range(NGEN):
+    print(f"-- Generation {g} --")
+    offspring = toolbox.select(pop, len(pop))
+    offspring = list(map(toolbox.clone, offspring))
+
+    # Crossover and mutation
+    for child1, child2 in zip(offspring[::2], offspring[1::2]):
+        if random.random() < CXPB:
+            toolbox.mate(child1, child2)
+            del child1.fitness.values
+            del child2.fitness.values
+
+    for mutant in offspring:
+        if random.random() < MUTPB:
+            toolbox.mutate(mutant)
+            del mutant.fitness.values
+
+    # Evaluate the individuals with an invalid fitness
+    invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
+    for ind in invalid_ind:
+        W = W_orig.copy()
+        specify_routes(W, ind)
+        ind.W = W
+    fitnesses = map(toolbox.evaluate, [ind.W for ind in invalid_ind])
+    for ind, fit in zip(invalid_ind, fitnesses):
+        ind.fitness.values = fit
+
+    # Print the best individual
+    best_ind = tools.selBest(pop, 1)[0]
+    print("")
+    print("Best individual: ", best_ind)
+    print("Fitness: ", best_ind.fitness.values[0])
+    display(best_ind.W.analyzer.basic_to_pandas())
+
+    # Update the population
+    pop[:] = offspring
+
+W_dso = best_ind.W.copy()
+
+##############################################################
+# Compare DUO and near-DSO
+
+print("DUO")
+W_duo.analyzer.macroscopic_fundamental_diagram()
+W_duo.analyzer.network_anim(file_name="out/grid_duo.gif", detailed=1, network_font_size=0, figsize=(6,6))
+
+print("near-DUO")
+W_dso.analyzer.macroscopic_fundamental_diagram()
+W_dso.analyzer.network_anim(file_name="out/grid_dso.gif", detailed=1, network_font_size=0, figsize=(6,6))
+
+print("Vehicle comparison")
+figure()
+subplot(111, aspect="equal")
+hist2d(
+    [veh.travel_time for veh in W_duo.VEHICLES.values()],
+    [veh.travel_time for veh in W_dso.VEHICLES.values()],
+    bins=20, range=[[0,1000],[0,1000]], cmap="Blues", cmin=1
+)
+colorbar().set_label("number of vehicles")
+plot([0,1000], [0,1000], "k--")
+xlabel("travel time of each vehicle in DUO (s)")
+ylabel("travel time of each vehicle in DSO (s)")
+show()

uxsim/Utilities/Utilities.py

@@ -2,5 +2,83 @@
 Submodule for general utilities.
 This (will) contains functions that are not essential for simulation but useful to specific analysis.
 """
+import networkx as nx
 
-pass
+def generate_grid_network(W, imax, jmax, **kwargs):
+    """
+    Generate a grid network with imax x jmax nodes.
+
+    Parameters
+    ----------
+    W : World
+        The world object to which the network will be added.
+    imax : int
+        The number of nodes in the x direction.
+    jmax : int
+        The number of nodes in the y direction. 
+    **kwargs : dict
+        Additional keyword arguments to be passed to the addLink function.
+    """
+
+    # Define the scenario
+    #deploy nodes as an imax x jmax grid
+    imax = 3
+    jmax = 3
+    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)
+
+    #create links between neighborhood nodes
+    links = {}
+    for i in range(imax):
+        for j in range(jmax):
+            if i != imax-1:
+                links[i,j,i+1,j] = W.addLink(f"l{(i,j,i+1,j)}", nodes[i,j], nodes[i+1,j], **kwargs)
+            if i != 0:
+                links[i,j,i-1,j] = W.addLink(f"l{(i,j,i-1,j)}", nodes[i,j], nodes[i-1,j], **kwargs)
+            if j != jmax-1:
+                links[i,j,i,j+1] = W.addLink(f"l{(i,j,i,j+1)}", nodes[i,j], nodes[i,j+1], **kwargs)
+            if j != 0:
+                links[i,j,i,j-1] = W.addLink(f"l{(i,j,i,j-1)}", nodes[i,j], nodes[i,j-1], **kwargs)
+
+    return nodes, links
+
+def enumerate_k_shortest_routes(W, source, target, k, cost_function=lambda l: l.length/l.u, print_stats=0):
+    """"
+    Enumerate the k shortest routes between two nodes in a network.
+
+    Parameters
+    ----------
+    W : World
+        The world object containing the network.   
+    source : str | Node
+        The source node.
+    target : str | Node
+        The target node.
+    k : int
+        The number of shortest routes to enumerate.
+    cost_function : function
+        A link cost function to compute shortest path. Default is the free-flow travel time.
+    print_stats : bool
+        Print the statistics of the paths.
+    """
+    G = nx.DiGraph()
+    for l in W.LINKS:
+        G.add_edge(l.start_node.name, l.end_node.name, weight=cost_function(l))
+
+    k_shortest_paths = list(nx.shortest_simple_paths(G, W.get_node(source).name, W.get_node(target).name, weight='weight'))[:k]
+    routes = []
+    for path in k_shortest_paths:
+        route = []
+        for n in path[:-1]:
+            for l in W.LINKS:
+                if l.start_node.name == n and l.end_node.name == path[path.index(n)+1]:
+                    route.append(l.name)
+        routes.append(route)
+
+        path_weight = sum([G[u][v]['weight'] for u, v in zip(path[:-1], path[1:])])
+        if print_stats:
+            print(f"Route: {route}, Cost: {path_weight}")
+
+    return routes

uxsim/analyzer.py

@@ -766,7 +766,7 @@ def show_simulation_progress(s):
             print(f"{s.W.TIME:>8.0f} s| {sum_vehs:>8.0f} vehs|  {avev:>4.1f} m/s| {time.time()-s.W.sim_start_time:8.2f} s", flush=True)
 
     @catch_exceptions_and_warn()
-    def network_anim(s, animation_speed_inverse=10, detailed=0, minwidth=0.5, maxwidth=12, left_handed=1, figsize=(6,6), node_size=2, network_font_size=20, timestep_skip=24):
+    def network_anim(s, animation_speed_inverse=10, detailed=0, minwidth=0.5, maxwidth=12, left_handed=1, figsize=(6,6), node_size=2, network_font_size=20, timestep_skip=24, file_name=None):
         """
         Generates an animation of the entire transportation network and its traffic states over time.
 
@@ -793,6 +793,8 @@ def network_anim(s, animation_speed_inverse=10, detailed=0, minwidth=0.5, maxwid
             The font size for the network labels in the animation. Default is 20.
         timestep_skip : int, optional
             How many timesteps are skipped per frame. Large value means coarse and lightweight animation. Default is 8.
+        file_name : str, optional
+            The name of the file to which the animation is saved. It overrides the defauld name. Default is None.
 
         Notes
         -----
@@ -812,12 +814,16 @@ def network_anim(s, animation_speed_inverse=10, detailed=0, minwidth=0.5, maxwid
                 else:
                     s.network_pillow(int(t), detailed=detailed, minwidth=minwidth, maxwidth=maxwidth, left_handed=left_handed, tmp_anim=1, figsize=figsize, node_size=node_size, network_font_size=network_font_size)
                 pics.append(Image.open(f"out{s.W.name}/tmp_anim_{t}.png"))
-        pics[0].save(f"out{s.W.name}/anim_network{detailed}.gif", save_all=True, append_images=pics[1:], optimize=False, duration=animation_speed_inverse*timestep_skip, loop=0)
+
+        fname = f"out{s.W.name}/anim_network{detailed}.gif"
+        if file_name != None:
+            fname = file_name
+        pics[0].save(fname, save_all=True, append_images=pics[1:], optimize=False, duration=animation_speed_inverse*timestep_skip, loop=0)
         for f in glob.glob(f"out{s.W.name}/tmp_anim_*.png"):
             os.remove(f)
 
     @catch_exceptions_and_warn()
-    def network_fancy(s, animation_speed_inverse=10, figsize=6, sample_ratio=0.3, interval=5, network_font_size=0, trace_length=3, speed_coef=2):
+    def network_fancy(s, animation_speed_inverse=10, figsize=6, sample_ratio=0.3, interval=5, network_font_size=0, trace_length=3, speed_coef=2, file_name=None):
         """
         Generates a visually appealing animation of vehicles' trajectories across the entire transportation network over time.
 
@@ -837,6 +843,8 @@ def network_fancy(s, animation_speed_inverse=10, figsize=6, sample_ratio=0.3, in
             The length of the vehicles' trajectory trails in the animation. Default is 3.
         speed_coef : int, optional
             A coefficient that adjusts the animation speed. Default is 2.
+        file_name : str, optional
+            The name of the file to which the animation is saved. It overrides the defauld name. Default is None.
 
         Notes
         -----
@@ -976,7 +984,10 @@ def flip(y):
             img.save(f"out{s.W.name}/tmp_anim_{t}.png")
             pics.append(Image.open(f"out{s.W.name}/tmp_anim_{t}.png"))
 
-        pics[0].save(f"out{s.W.name}/anim_network_fancy.gif", save_all=True, append_images=pics[1:], optimize=False, duration=animation_speed_inverse*speed_coef, loop=0)
+        fname = f"out{s.W.name}/anim_network_fancy.gif"
+        if file_name != None:
+            fname = file_name
+        pics[0].save(fname, save_all=True, append_images=pics[1:], optimize=False, duration=animation_speed_inverse*speed_coef, loop=0)
 
         for f in glob.glob(f"out{s.W.name}/tmp_anim_*.png"):
             os.remove(f)

uxsim/utils.py

@@ -132,6 +132,19 @@ def print_columns(*lists):
                 print(" " * 10, end=" ")  # Adjust spacing to match the above width
         print()  # Newline after each row
 
+def display_image_in_notebook(image_path):
+    """
+    Display an image in Jupyter Notebook.
+
+    Parameters
+    ----------
+    image_path : str
+        The path to the image file to display.
+    """
+    from IPython.display import display, Image
+    with open(image_path, "rb") as f:
+        display(Image(data=f.read(), format='png'))
+
 class LoggingWarning(UserWarning):
     """
     This warns that when vehicle_logging_timestep_interval is not 1 but called vehicle logging-related functions.