Route choice docs and adjusting

aequilibrae/paths/route_choice.py

@@ -222,7 +222,7 @@ def execute_single(self, origin: int, destination: int, perform_assignment: bool
             **self.parameters,
         )
 
-    def execute(self, perform_assignment: bool = False) -> None:
+    def execute(self, perform_assignment: bool = True) -> None:
         """
         Generate route choice sets between the previously supplied nodes, potentially performing an assignment.
 

docs/source/conf.py

@@ -77,7 +77,7 @@
             "examples/trip_distribution",
             "examples/visualization",
             "examples/aequilibrae_without_a_model",
-            "examples/full_workflows",
+            "examples/assignment_workflows",
             "examples/other_applications",
         ]
     ),

docs/source/examples/trip_distribution/plot_route_choice.py → docs/source/examples/assignment_workflows/plot_route_choice.py

@@ -7,12 +7,16 @@
 Serena Metropolitan Area in Chile.
 
 """
+# %%
 
 # Imports
 from uuid import uuid4
 from tempfile import gettempdir
 from os.path import join
 from aequilibrae.utils.create_example import create_example
+# sphinx_gallery_thumbnail_path = 'images/plot_route_choice_assignment.png'
+
+# %%
 
 # We create the example project inside our temp folder
 fldr = join(gettempdir(), uuid4().hex)
@@ -37,6 +41,19 @@
 # %%
 import numpy as np
 
+# %%
+# Model parameters
+# ~~~~~~~~~~~~~~~~
+# We'll set the parameters for our route choice model. These are the parameters that will be used to calculate the
+# utility of each path. In our example, the utility is equal to *theta* * distance
+# And the path overlap factor (PSL) is equal to *beta*.
+
+# Distance factor
+theta=0.011
+
+# PSL parameter
+beta = 1.1
+
 # %%
 # Let's build all graphs
 project.network.build_graphs()
@@ -50,12 +67,18 @@
 # we also see what graphs are available
 project.network.graphs.keys()
 
-# let's say we want to minimize the distance
-graph.set_graph("distance")
+od_pairs_of_interest = [(71645, 79385), (77011, 74089)]
+nodes_of_interest = (71645, 74089, 77011, 79385)
 
-# But let's say we only want a skim matrix for nodes 28-40, and 49-60 (inclusive), these happen to be a selection of
-# western centroids.
-graph.prepare_graph(np.array(list(range(28, 41)) + list(range(49, 91))))
+# let's say that utility is just a function of distance
+# So we build our *utility* field as the distance times theta
+graph.network = graph.network.assign(utility=graph.network.distance * theta)
+
+# Prepare the graph with all nodes of interest as centroids
+graph.prepare_graph(np.array(nodes_of_interest))
+
+# And set the cost of the graph the as the utility field just created
+graph.set_graph("utility")
 
 # %%
 # Mock demand matrix
@@ -68,8 +91,8 @@
 mat = AequilibraeMatrix()
 mat.create_empty(zones=graph.num_zones, matrix_names=names_list, memory_only=True)
 mat.index = graph.centroids[:]
-mat.matrices[:, :, 0] = np.full((graph.num_zones, graph.num_zones), 1.0)
-mat.matrices[:, :, 1] = np.full((graph.num_zones, graph.num_zones), 5.0)
+mat.matrices[:, :, 0] = np.full((graph.num_zones, graph.num_zones), 10.0)
+mat.matrices[:, :, 1] = np.full((graph.num_zones, graph.num_zones), 50.0)
 mat.computational_view()
 
 # %%
@@ -109,8 +132,8 @@
 #
 # It is highly recommended to set either `max_routes` or `max_depth` to prevent runaway results.
 
-# rc.set_choice_set_generation("link-penalisation", max_routes=5, penalty=1.1)
-rc.set_choice_set_generation("bfsle", max_routes=5, beta=1.1, theta=1.1)
+# rc.set_choice_set_generation("link-penalisation", max_routes=5, penalty=1.02)
+rc.set_choice_set_generation("bfsle", max_routes=5, beta=1.1)
 
 # %%
 # All parameters are optional, the defaults are:
@@ -119,14 +142,53 @@
 # %%
 # We can now perform a computation for single OD pair if we'd like. Here we do one between the first and last centroid
 # as well an an assignment.
-results = rc.execute_single(28, 90, perform_assignment=True)
+results = rc.execute_single(77011, 74089, perform_assignment=True)
 print(results[0])
 
 # %%
 # Because we asked it to also perform an assignment we can access the various results from that
 # The default return is a Pyarrow Table but Pandas is nicer for viewing.
 rc.get_results().to_pandas()
 
+
+# %%
+# let's define a function to plot assignment results
+
+def plot_results(link_loads):
+    import folium
+    import geopandas as gpd
+
+
+    link_loads =link_loads[["link_id", "demand_tot"]]
+    max_load = link_loads["demand_tot"].max()
+    links = gpd.GeoDataFrame(project.network.links.data, crs=4326)
+    links = links.merge(link_loads, on="link_id")
+    links = links[links["demand_tot"] > 0]
+
+    loads_lyr = folium.FeatureGroup("link_loads")
+
+    # Maximum thickness we would like is probably a 10, so let's make sure we don't go over that
+    factor = 10 / max_load
+
+    # Let's create the layers
+    for _, rec in links.iterrows():
+        points = rec.geometry.wkt.replace("LINESTRING ", "").replace("(", "").replace(")", "").split(", ")
+        points = "[[" + "],[".join([p.replace(" ", ", ") for p in points]) + "]]"
+        # we need to take from x/y to lat/long
+        points = [[x[1], x[0]] for x in eval(points)]
+        _ = folium.vector_layers.PolyLine(points, color="red", weight=factor * rec.demand_tot).add_to(loads_lyr)
+
+    long, lat = project.conn.execute("select avg(xmin), avg(ymin) from idx_links_geometry").fetchone()
+
+    map_osm = folium.Map(location=[lat, long], tiles="Cartodb Positron", zoom_start=12)
+    loads_lyr.add_to(map_osm)
+    folium.LayerControl().add_to(map_osm)
+    return map_osm
+
+# %%
+plot_results(rc.get_load_results()[0])
+
+
 # %%
 # To perform a batch operation we need to prepare the object first. We can either provide a list of tuple of the OD
 # pairs we'd like to use, or we can provided a 1D list and the generation will be run on all permutations.
@@ -141,6 +203,10 @@
 # Since we provided a matrix initially we can also perform link loading based on our assignment results.
 rc.get_load_results()
 
+# %% we can plot these as well
+plot_results(rc.get_load_results()[0])
+
+
 # %%
 # Select link analysis
 # ~~~~~~~~~~~~~~~~~~

docs/source/examples/assignment_workflows/plot_route_choice_set.py

@@ -0,0 +1,166 @@
+""".. _example_usage_route_choice_generation:
+
+Route Choice set generation
+===========================
+
+In this example, we show how to generate route choice sets for estimation of route choice models, using a
+a city in La Serena Metropolitan Area in Chile.
+
+"""
+
+# %%
+
+# Imports
+from uuid import uuid4
+from tempfile import gettempdir
+from os.path import join
+import numpy as np
+from aequilibrae.utils.create_example import create_example
+# sphinx_gallery_thumbnail_path = 'images/plot_route_choice_set.png'
+
+# %%
+# We create the example project inside our temp folder
+fldr = join(gettempdir(), uuid4().hex)
+
+project = create_example(fldr, "coquimbo")
+
+# %%
+# Choice set generation
+# ---------------------
+
+
+# %%
+
+od_pairs_of_interest = [(71645, 79385), (77011, 74089)]
+nodes_of_interest = (71645, 74089, 77011, 79385)
+
+# %%
+# Let's build all graphs
+project.network.build_graphs()
+# We get warnings that several fields in the project are filled with NaNs.
+# This is true, but we won't use those fields.
+
+
+# %%
+# We grab the graph for cars
+graph = project.network.graphs["c"]
+
+# we also see what graphs are available
+project.network.graphs.keys()
+
+graph.set_graph("distance")
+
+# We set the nodes of interest as centroids to make sure they are not simplified away when we create the network
+graph.prepare_graph(np.array(nodes_of_interest))
+
+# We allow flows through "centroid connectors" because our centroids are not really centroids
+# If we have actual centroid connectors in the network (and more than one per centroid) , then we
+# should remove them from the graph
+graph.set_blocked_centroid_flows(False)
+
+# %%
+# Route Choice class
+# ~~~~~~~~~~~~~~~~~~
+# Here we'll construct and use the Route Choice Set class to generate our route sets
+from aequilibrae.paths.route_choice_set import RouteChoiceSet
+
+# %%
+# This object construct might take a minute depending on the size of the graph due to the construction of the compressed
+# link to network link mapping that's required.  This is a one time operation per graph and is cached. We need to
+# supply a Graph and optionally a AequilibraeMatrix, if the matrix is not provided link loading cannot be preformed.
+rc = RouteChoiceSet(graph)
+
+# %%
+
+# Here we'll set the parameters of our set generation. There are two algorithms available: Link penalisation, and BFSLE
+# based on the paper
+# "Route choice sets for very high-resolution data" by Nadine Rieser-Schüssler, Michael Balmer & Kay W. Axhausen (2013).
+# https://doi.org/10.1080/18128602.2012.671383
+#
+# Our BFSLE implementation has been extended to allow applying link penalisation as well. Every
+# link in all routes found at a depth are penalised with the `penalty` factor for the next depth. So at a depth of 0 no
+# links are penalised nor removed. At depth 1, all links found at depth 0 are penalised, then the links marked for
+# removal are removed. All links in the routes found at depth 1 are then penalised for the next depth. The penalisation
+# compounds. Pass set `penalty=1.0` to disable.
+#
+# It is highly recommended to set either `max_routes` or `max_depth` to prevent runaway results.
+
+# rc.set_choice_set_generation("link-penalisation", max_routes=5, penalty=1.02)
+# The 5% penalty (1.05) is likely a little too large, but it create routes that are distinct enough to make this simple
+# example more interesting
+rc.batched(od_pairs_of_interest, max_routes=5, cores=10, bfsle=False, penalty=1.05, path_size_logit=True)
+choice_set = rc.get_results().to_pandas()
+
+# %%
+# Plotting choice sets
+# --------------------
+
+
+# %%
+# Now we will plot the paths we just created for the second OD pair
+import folium
+import geopandas as gpd
+
+# %%
+
+# Let's create a separate for each route so we can visualize one at a time
+rlyr1 = folium.FeatureGroup("route 1")
+rlyr2 = folium.FeatureGroup("route 2")
+rlyr3 = folium.FeatureGroup("route 3")
+rlyr4 = folium.FeatureGroup("route 4")
+rlyr5 = folium.FeatureGroup("route 5")
+od_lyr = folium.FeatureGroup("Origin and Destination")
+layers = [rlyr1, rlyr2, rlyr3, rlyr4, rlyr5]
+
+# %%
+
+# We get the data we will use for the plot: Links, Nodes and the route choice set
+links = gpd.GeoDataFrame(project.network.links.data, crs=4326)
+nodes = gpd.GeoDataFrame(project.network.nodes.data, crs=4326)
+
+plot_routes = choice_set[(choice_set["origin id"] == 77011)]["route set"].values
+
+# Let's create the layers
+colors = ["red", "blue", "green", "purple", "orange"]
+for i, route in enumerate(plot_routes):
+    rt = links[links.link_id.isin(route)]
+    routes_layer = layers[i]
+    for wkt in rt.geometry.to_wkt().values:
+        points = wkt.replace("LINESTRING ", "").replace("(", "").replace(")", "").split(", ")
+        points = "[[" + "],[".join([p.replace(" ", ", ") for p in points]) + "]]"
+        # we need to take from x/y to lat/long
+        points = [[x[1], x[0]] for x in eval(points)]
+
+        _ = folium.vector_layers.PolyLine(points, color=colors[i], weight=4).add_to(routes_layer)
+
+# Creates the points for both origin and destination
+for i, row in nodes[nodes.node_id.isin((77011, 74089))].iterrows():
+    point = (row.geometry.y, row.geometry.x)
+
+    _ = folium.vector_layers.CircleMarker(
+        point,
+        popup=f"<b>link_id: {row.node_id}</b>",
+        color="red",
+        radius=5,
+        fill=True,
+        fillColor="red",
+        fillOpacity=1.0,
+    ).add_to(od_lyr)
+
+# %%
+# It is worthwhile to notice that using distance as the cost function, the routes are not the fastest ones as the
+# freeway does not get used
+
+# %%
+# Create the map and center it in the correct place
+long, lat = project.conn.execute("select avg(xmin), avg(ymin) from idx_links_geometry").fetchone()
+
+map_osm = folium.Map(location=[lat, long], tiles="Cartodb Positron", zoom_start=12)
+for routes_layer in layers:
+    routes_layer.add_to(map_osm)
+od_lyr.add_to(map_osm)
+folium.LayerControl().add_to(map_osm)
+map_osm
+
+# %%
+project.close()

docs/source/examples/assignment_workflows/readme.rst

@@ -0,0 +1,2 @@
+Assignment Workflows
+--------------------