Skip to main content
`Ctrl`+`K`

Note

Go to the end to download the full example code.

# Route Choice

In this example, we show how to perform route choice set generation using BFSLE and Link penalisation, for a city in La Serena Metropolitan Area in Chile.

References

- Route Choice

See also

Several functions, methods, classes and modules are used in this example:

- `aequilibrae.paths.Graph()`

- `aequilibrae.paths.RouteChoice()`

- `aequilibrae.matrix.AequilibraeMatrix()`

    # Imports
    from uuid import uuid4
    from tempfile import gettempdir
    from os.path import join
    from aequilibrae.utils.create_example import create_example

    # We create the example project inside our temp folder
    fldr = join(gettempdir(), uuid4().hex)

    project = create_example(fldr, "coquimbo")

    import logging
    import sys

    # When the project opens, we can tell the logger to direct all messages to the terminal as well
    logger = project.logger
    stdout_handler = logging.StreamHandler(sys.stdout)
    formatter = logging.Formatter("%(asctime)s;%(levelname)s ; %(message)s")
    stdout_handler.setFormatter(formatter)
    logger.addHandler(stdout_handler)

## Model parameters

    import numpy as np

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 \(distance * theta\), and the path overlap factor (PSL) is equal to \(beta\).

    # Distance factor
    theta = 0.00011

    # PSL parameter
    beta = 1.1

Let’s select a set of nodes of interest

    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.

    2024-11-04 00:39:33,040;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:39:33,116;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:39:33,208;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations
    2024-11-04 00:39:33,300;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations

We also see what graphs are available

    project.network.graphs.keys()

    dict_keys(['b', 'c', 't', 'w'])

We grab the graph for cars

    graph = project.network.graphs["c"]

    # Let's say that utility is just a function of distance, so we build our 'utility' field as distance * 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")

    # 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)

    2024-11-04 00:39:33,366;WARNING ; Field(s) speed, travel_time, capacity, osm_id, lanes has(ve) at least one NaN value. Check your computations

## Mock demand matrix

We’ll create a mock demand matrix with demand 1 for every zone and prepare it for computation.

    from aequilibrae.matrix import AequilibraeMatrix

    names_list = ["demand", "5x demand"]

    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), 10.0)
    mat.matrices[:, :, 1] = np.full((graph.num_zones, graph.num_zones), 50.0)
    mat.computational_view()

## Create plot function

Before dive into the Route Choice class, let’s define a function to plot assignment results.

    import folium

    def plot_results(link_loads):

        link_loads = link_loads[link_loads.tot > 0]
        max_load = link_loads["tot"].max()
        links = project.network.links.data
        loaded_links = links.merge(link_loads, on="link_id", how="inner")

        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 loaded_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,
                tooltip=f"link_id: {rec.link_id}, Flow: {rec.tot:.3f}",
                color="red",
                weight=factor * rec.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

## Route Choice class

Here we’ll construct and use the Route Choice class to generate our route sets

    from aequilibrae.paths import RouteChoice

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.

    rc = RouteChoice(graph)

    # Let's check the default parameters for the Route Choice class
    print(rc.default_parameters)

    {'generic': {'seed': 0, 'max_routes': 0, 'max_depth': 0, 'max_misses': 100, 'penalty': 1.01, 'cutoff_prob': 0.0, 'beta': 1.0, 'store_results': True}, 'link-penalisation': {}, 'bfsle': {'penalty': 1.0}}

Let’s add the demand. If it’s not provided, link loading cannot be preformed.

    rc.add_demand(mat)

It is highly recommended to set either `max_routes` or `max_depth` to prevent runaway results.

    rc.set_choice_set_generation("bfsle", max_routes=5)

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 as an assignment.

    results = rc.execute_single(77011, 74089, demand=1.0)
    print(results[0])

    (24222, 30332, 30333, 10435, 30068, 30069, 14198, 14199, 31161, 30928, 30929, 30930, 30931, 24172, 30878, 30879, 30880, 30881, 30882, 30883, 30884, 30885, 30886, 30887, 30888, 30889, 30890, 30891, 5179, 5180, 5181, 5182, 26463, 26462, 26461, 26460, 26459, 26458, 26457, 26456, 26480, 3341, 3342, 3339, 9509, 9510, 9511, 9512, 18487, 14972, 14973, 32692, 32693, 32694, 2300, 2301, 33715, 19978, 19979, 19977, 19976, 19975, 19974, 19973, 19972, 19971, 19970, 22082, 22080, 5351, 5352, 2280, 2281, 2282, 575, 576, 577, 578, 579, 536, 537, 538, 539, 540, 541, 15406, 15407, 15408, 553, 552, 633, 634, 635, 630, 631, 632, 623, 624, 625, 626, 471, 5363, 34169, 34170, 34171, 34785, 6466, 6465, 29938, 29939, 29940, 29941, 1446, 1447, 1448, 1449, 1450, 939, 940, 941, 9840, 9841, 26314, 26313, 26312, 26311, 26310, 26309, 26308, 26307, 26306, 26305, 26304, 26303, 26302, 26301, 26300, 34079, 34147, 29962, 26422, 26421, 26420, 765, 764, 763, 762, 761, 760, 736, 10973, 10974, 10975, 725, 10972, 727, 728, 26424, 733, 734, 29899, 20970, 20969, 20968, 20967, 20966, 20965, 20964, 20963, 20962, 9584, 9583, 20981, 21398, 20982, 34208, 35, 36, 59, 60, 61, 22363, 22364, 22365, 22366, 22367, 28958, 28959, 28960, 28961, 28962, 28805, 28806, 28807, 28808, 28809, 28810, 28827, 28828, 28829, 28830, 28874)

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.

    res = rc.get_results().to_pandas()
    res.head()

| origin id | destination id | route set | cost | mask | path overlap | probability
--- | --- | --- | --- | --- | --- | --- | ---
0 | 77011 | 74089 | [24222, 30332, 30333, 10435, 30068, 30069, 141... | 1.907420 | True | 0.386507 | 0.231745
1 | 77011 | 74089 | [24222, 30332, 30333, 10435, 30068, 30069, 141... | 1.890617 | True | 0.287132 | 0.175078
2 | 77011 | 74089 | [24223, 31635, 31636, 31637, 31638, 31639, 316... | 1.896528 | True | 0.304541 | 0.184599
3 | 77011 | 74089 | [24222, 30332, 30333, 10435, 30068, 30069, 141... | 1.909910 | True | 0.430542 | 0.257506
4 | 77011 | 74089 | [24222, 30332, 30333, 10435, 30068, 30069, 141... | 1.886882 | True | 0.246836 | 0.151071

    plot_results(rc.get_load_results()["demand"])

Make this Notebook Trusted to load map: File -> Trust Notebook

## Batch operations

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.

    rc.prepare()

Now we can perform a batch computation with an assignment

    rc.execute(perform_assignment=True)
    res = rc.get_results().to_pandas()
    res.head()

| origin id | destination id | route set | cost | mask | path overlap | probability
--- | --- | --- | --- | --- | --- | --- | ---
0 | 71645 | 74089 | [19550, 19549, 19548, 19547, 19546, 19545, 209... | 0.607017 | True | 0.384514 | 0.204478
1 | 71645 | 74089 | [19550, 19549, 19548, 19547, 19546, 19545, 195... | 0.607425 | True | 0.273951 | 0.145623
2 | 71645 | 74089 | [19550, 19549, 19548, 19547, 19546, 19545, 209... | 0.622079 | True | 0.544362 | 0.285155
3 | 71645 | 74089 | [19550, 19549, 19548, 19547, 19546, 19545, 209... | 0.611944 | True | 0.443492 | 0.234682
4 | 71645 | 74089 | [19550, 19549, 19548, 19547, 19546, 19545, 209... | 0.604384 | True | 0.243935 | 0.130063

Since we provided a matrix initially we can also perform link loading based on our assignment results.

    rc.get_load_results()

| demand | 5x demand | demand | 5x demand
--- | --- | --- | --- | ---
| ab | ba | ab | ba | tot | tot
link_id |  |  |  |  |  |
1 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
2 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
3 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
12 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
13 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
... | ... | ... | ... | ... | ... | ...
34938 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
34939 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
34940 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
34941 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
34942 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0

19983 rows × 6 columns

We can plot these as well

    plot_results(rc.get_load_results()["demand"])

Make this Notebook Trusted to load map: File -> Trust Notebook

## Select link analysis

We can also enable select link analysis by providing the links and the directions that we are interested in. Here we set the select link to trigger when (7369, 1) and (20983, 1) is utilised in “sl1” and “sl2” when (7369, 1) is utilised.

    rc.set_select_links({"sl1": [[(7369, 1), (20983, 1)]], "sl2": [[(7369, 1)]]})
    rc.execute(perform_assignment=True)

We can get then the results in a Pandas DataFrame for both the network.

    sl = rc.get_select_link_loading_results()
    sl

| demand | 5x demand | demand | 5x demand | demand | 5x demand | demand | 5x demand
--- | --- | --- | --- | --- | --- | --- | --- | ---
| sl1 | sl1 | sl2 | sl2 | sl1 | sl1 | sl2 | sl2
| ab | ba | ab | ba | ab | ba | ab | ba | tot | tot | tot | tot
link_id |  |  |  |  |  |  |  |  |  |  |  |
1 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
2 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
3 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
12 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
13 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ...
34938 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
34939 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
34940 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
34941 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
34942 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0

19983 rows × 12 columns

We can also access the OD matrices for this link loading. These matrices are sparse and can be converted to SciPy sparse matrices for ease of use. They’re stored in a dictionary where the key is the matrix name concatenated with the select link set name via an underscore.

    rc.get_select_link_od_matrix_results()

    {'sl1': {'demand': <aequilibrae.matrix.sparse_matrix.COO object at 0x7eff5bea2680>, '5x demand': <aequilibrae.matrix.sparse_matrix.COO object at 0x7eff5bea3b20>}, 'sl2': {'demand': <aequilibrae.matrix.sparse_matrix.COO object at 0x7eff5bea02e0>, '5x demand': <aequilibrae.matrix.sparse_matrix.COO object at 0x7eff5bea2500>}}

    od_matrix = rc.get_select_link_od_matrix_results()["sl1"]["demand"]
    od_matrix.to_scipy().toarray()

    array([[0.        , 0.        , 0.        , 3.04610785],
           [0.        , 0.        , 0.        , 0.        ],
           [0.        , 0.        , 0.        , 0.        ],
           [0.        , 0.        , 0.        , 0.        ]])

    project.close()

    INFO:aequilibrae:Closed project on /tmp/4797712ab95d43959090349949f8818e

Total running time of the script: (0 minutes 4.103 seconds)

Gallery generated by Sphinx-Gallery