Extra material: Shortest path problem in real life#
import os
import sys
import subprocess
try:
import google.colab
from google.colab import files
IN_COLAB = True
except ImportError:
IN_COLAB = False
if IN_COLAB:
subprocess.check_call([
sys.executable, "-m", "pip", "install", "-q",
"pandana==0.8", "osmnx", "geopy", "folium", "matplotlib",
"scikit-learn", "psutil", "mapclassify", "threadpoolctl",
])
# Set before importing Pandana, so OpenMP initializes with one worker.
# Deliberately compare both routing implementations without parallel speedup.
PANDANA_THREADS = 1
os.environ["OMP_NUM_THREADS"] = str(PANDANA_THREADS)
os.environ["OMP_DYNAMIC"] = "FALSE"
import gc
import zipfile
from numbers import Number
import hashlib
import json
import platform
from datetime import datetime, timezone
from importlib.metadata import version
from pathlib import Path
from time import perf_counter
import folium
import geopy
from geopy.extra.rate_limiter import RateLimiter
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import osmnx as ox
import pandas as pd
import pandana
import psutil
from threadpoolctl import threadpool_info
from IPython.display import display, FileLink
from packaging.version import Version
assert Version(version("pandana")) >= Version("0.8"), "Please upgrade Pandana."
pd.options.display.float_format = "{:.3f}".format
# The default runs without Photon or Overpass requests.
RUN_LIVE = False
RUN_FULL_BENCHMARK = False
LIVE_TIMEOUT_S = 180
evidence = Path("data") / "shortest-path-live"
evidence.mkdir(parents=True, exist_ok=True)
openmp_pools = [p for p in threadpool_info() if p["user_api"] == "openmp"]
assert openmp_pools and all(p["num_threads"] == 1 for p in openmp_pools)
def style_table(table):
"""Align text left and numbers right, including row and column headings."""
def alignment(value):
return "text-align: right" if isinstance(value, Number) else "text-align: left"
styles = [
{"selector": f"th.col_heading.col{i}",
"props": [("text-align", "right" if pd.api.types.is_numeric_dtype(table[column])
else "left")]}
for i, column in enumerate(table.columns)
]
return (table.style.map(alignment)
.map_index(alignment, axis=0)
.set_table_styles(styles))
Introduction#
Google Maps brought the world to our screens, including geocoding and routing for several modes of transport. Much of its use is interactive, but data and analytics applications also need programmatic access to these services. Here we use open tools to geocode places and find walking routes in Amsterdam.
The main question is how to compute many shortest paths efficiently. We begin with individual routes, then compare a distance matrix computed row by row with one-to-many Dijkstra searches against a single batched Pandana call. Pandana uses contraction hierarchies to preprocess the network for repeated queries.
For local execution install pandana>=0.8, osmnx, geopy, folium,
matplotlib, scikit-learn, mapclassify, psutil, and threadpoolctl. The recorded experiment was tested with Pandana 0.8; Colab installation is handled above.
The default uses the bundled Amsterdam graph and geocoding responses. Interactive
maps still need internet access for their JavaScript and background tiles. If the notebook is opened alone on
Colab, it downloads the fixed input archive from GitHub with a checksum check.
Running and saving in Google Colab#
Open this notebook in Colab with a fresh CPU runtime, then choose Runtime >
Run all. The first cell installs Pandana 0.8 and the required packages. The
default uses the recorded graph and the six recorded place responses together.
Set RUN_LIVE = True to refresh both from Photon and OpenStreetMap. Individual
requests have timeouts, and the complete refresh runs in a subprocess with a
180-second limit, including retries and server waiting. A failed refresh falls
back to the complete recorded pair and prints that choice; live and recorded
inputs are never mixed. Live results can change with map updates. The benchmark records the graph checksum, sampled nodes,
installed versions and the actual runtime CPU and memory. The standard NetworkX
Dijkstra implementation used here is single-threaded, and we process origins
sequentially. We deliberately restrict Pandana to one thread to showcase the
difference between repeated Dijkstra searches and preprocessed
contraction-hierarchy queries without a parallel speedup advantage.
NetworkX uses Python and Pandana’s routing engine uses C++, so language and
implementation overhead also contribute; this is not a language-controlled test.
This describes the standard NetworkX routine used in this notebook, not a
limitation of every NetworkX backend.
PANDANA_THREADS sets OMP_NUM_THREADS before Pandana is imported, controlling
the OpenMP workers used for construction and batched queries. Here it is fixed at
1 throughout. OMP_DYNAMIC=FALSE disables dynamic adjustment, and the notebook
checks the loaded OpenMP runtime. Restart the runtime before rerunning the setup
cell if Pandana was already imported with another thread setting.
This single-thread setting is a deliberate choice for the comparison, not a recommendation for deployment. In your own application, allowing Pandana to use multiple threads may further reduce construction and batched-query times, depending on the workload and available CPU cores. Practical performance may therefore improve beyond the single-thread results shown here; benchmark the thread count on your target machine rather than assuming proportional speedup.
The final export cell downloads an archive of this session’s timing records and environment. Also use File > Download > Download .ipynb to keep the notebook with its rendered outputs; leave Omit code cell output when saving disabled. Runtime files alone are temporary. See the Colab FAQ. The national results remain explicitly labelled as previously measured Windows results; running Colab does not turn them into Colab measurements.
The default RUN_FULL_BENCHMARK = False measures matrices for 10, 25 and 50
sampled road nodes, with the same warm-up and three repetitions per method.
Set it to True for the original full sequence through 1,000 nodes. The full
NetworkX benchmark can take tens of minutes. Neither setting reruns the
historical national experiment. Generated session files are kept in an ignored
data/shortest-path-live/ folder.
Some background information and history#
Geoff Boeing has made urban data analytics more accessible, particularly the analysis of street networks. His research is accompanied by usable demonstrations with OSMnx.
Professor Peter Sanders and his collaborators’ route-planning research helped make exact shortest-path queries practical on very large road networks. Two contributors worth highlighting are:
Dominik Schultes, whose work with Sanders on highway hierarchies helped exploit the hierarchical structure of road networks. The subsequent work on contraction hierarchies, with Robert Geisberger and Daniel Delling, made preprocessing simpler while keeping queries fast. Exact shortest-path algorithms already existed; the advance was their practical scalability for repeated queries on large road networks. For historical context, see the The Shortest Path Problem (DIMACS Series, vol. 74, 2009) and the Scientific American 50 feature discussing Schultes and transit-node routing.
Dennis Luxen, a creator of the Open Source Routing Machine (OSRM), which makes scalable road routing available as open-source software.
Pandana originated in research led by Paul Waddell at the University of California, Berkeley, on land use, transportation, and accessibility. Waddell conceived the application with Fletcher Foti, his doctoral student, and Dennis Luxen, and provided research direction and funding. Foti led the original software implementation and its pandas-oriented interface; Luxen contributed routing expertise and contraction-hierarchy code. Their 2012 paper describes how the framework connects parcels and complete street networks for accessibility analysis from pedestrian to metropolitan scales.
Pandana subsequently became part of the Urban Data Science Toolkit (UDST), with contributions from Matt Davis, Federico Fernandez, Sam Maurer, and the wider community.
Pandana 0.8, released in August 2026, reflects a collective effort to support the modern Python ecosystem. Paul Waddell helped coordinate the renewed maintenance work, and Sam Maurer helped bring the release to completion (Paul Waddell, personal communication, August 2026). Among the wider community of contributors, Joaquim Gromicho contributed compatibility updates supporting Python 3.14, NumPy 2, and pandas 3. This contribution was part of the collective maintenance and release effort.
The key idea of contraction hierarchies is to remove nodes in a preprocessing order, adding shortcuts where needed to preserve shortest-path distances. Queries can then search the hierarchy rather than explore the full network. Preprocessing takes time and memory, but its cost can be shared across many queries. The experiment at the end examines this tradeoff.
Geocoding and map visualization#
The world is mapped with the geographic coordinate system but we have difficulties remembering latitudes and longitudes. We learn and remember the world better from addresses.
# The recorded graph and six geocoding responses always travel together.
import io
import tempfile
from urllib.request import urlopen
RECORDED_ARCHIVE_SHA256 = "b9bd52209bbacff3c956de5d08431f8dbc2155c4b30fe4862a1eaefa7c30fc3a"
RECORDED_URL = (
"https://raw.githubusercontent.com/mobook/MO-book/"
"e745704ca19c0eef476beb32bb24dfba80244043/"
"notebooks/04/data/road-network-recorded-results.zip"
)
def recorded_inputs(output_dir):
local = next((p for p in (
Path("data/road-network-recorded-results.zip"),
Path("notebooks/04/data/road-network-recorded-results.zip"),
) if p.is_file()), None)
if local is None:
with urlopen(RECORDED_URL, timeout=30) as response:
content = response.read()
else:
content = local.read_bytes()
if hashlib.sha256(content).hexdigest() != RECORDED_ARCHIVE_SHA256:
raise ValueError("Recorded archive checksum mismatch; restore the bundled snapshot.")
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
with zipfile.ZipFile(io.BytesIO(content)) as bundle:
graph = bundle.read("amsterdam-walk.graphml.gz")
geocoding = bundle.read("geocoding.json")
old_manifest = json.loads(bundle.read("manifest.json"))
graph_hash = hashlib.sha256(graph).hexdigest()
if graph_hash != old_manifest["amsterdam_environment"]["graph_sha256"]:
raise ValueError("Recorded graph checksum mismatch.")
(output_dir / "amsterdam-walk.graphml.gz").write_bytes(graph)
(output_dir / "geocoding.json").write_bytes(geocoding)
return output_dir, {
"mode": "recorded", "archive_sha256": RECORDED_ARCHIVE_SHA256,
"graph_sha256": graph_hash,
"geocoding_retrieved_utc": json.loads(geocoding)["retrieved_utc"],
"graph_recorded_run_utc": old_manifest["amsterdam_environment"]["run_utc"],
}
LIVE_PROGRAM = r'''
import json, sys
from pathlib import Path
from datetime import datetime, timezone
import geopy
from geopy.extra.rate_limiter import RateLimiter
import osmnx as ox
folder = Path(sys.argv[1])
names = ["Centraal Station", "VU Hoofdgebouw", "Amsterdam Business School",
"University of Amsterdam", "Johan Cruijff ArenA", "Ziggo Dome"]
locator = geopy.Photon(user_agent="MO-book-shortest-path-example", timeout=15)
geocode = RateLimiter(locator.geocode, min_delay_seconds=1, max_retries=1,
error_wait_seconds=3, swallow_exceptions=False)
locations = [geocode(name + ", Amsterdam, NL") for name in names]
if any(location is None for location in locations):
raise ValueError("At least one live place lookup returned no result.")
record = {"retrieved_utc": datetime.now(timezone.utc).isoformat(),
"provider": "Photon", "results": [loc.raw for loc in locations]}
ox.settings.log_console = False
ox.settings.requests_timeout = 30
ox.settings.use_cache = True
ox.settings.cache_folder = folder / ".cache"
graph = ox.graph_from_place("Amsterdam, NL", network_type="walk")
ox.save_graphml(graph, folder / "amsterdam-walk.graphml.gz")
(folder / "geocoding.json").write_text(json.dumps(record, indent=2), encoding="utf8")
'''
def acquire_inputs(output_dir, use_live=False, timeout_s=180):
"""Acquire a complete pair; a failed live refresh falls back as a whole."""
if not np.isfinite(timeout_s) or timeout_s <= 0:
raise ValueError("The live acquisition time limit must be positive and finite.")
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
failure = None
if use_live:
# A subprocess bounds OSMnx's retries and server-slot waiting as well as HTTP calls.
live_dir = Path(tempfile.mkdtemp(prefix="live-", dir=output_dir))
try:
subprocess.run([sys.executable, "-c", LIVE_PROGRAM, str(live_dir.resolve())],
timeout=timeout_s, check=True, capture_output=True, text=True)
graph_file = live_dir / "amsterdam-walk.graphml.gz"
record = json.loads((live_dir / "geocoding.json").read_text(encoding="utf8"))
if len(record["results"]) != 6:
raise ValueError("Expected all six live place responses.")
return live_dir, {
"mode": "live", "geocoding_retrieved_utc": record["retrieved_utc"],
"graph_sha256": hashlib.sha256(graph_file.read_bytes()).hexdigest(),
"acquired_utc": datetime.now(timezone.utc).isoformat(),
}
except subprocess.TimeoutExpired:
failure = f"Live acquisition exceeded {timeout_s:g} seconds."
except (subprocess.CalledProcessError, OSError, ValueError, KeyError) as error:
failure = f"Live acquisition failed ({type(error).__name__})."
print(f"{failure} Using the complete recorded graph/geocoding pair.")
folder, provenance = recorded_inputs(output_dir / "recorded-inputs")
if failure is not None:
provenance["live_failure"] = failure
return folder, provenance
inputs_dir, input_provenance = acquire_inputs(evidence, RUN_LIVE, LIVE_TIMEOUT_S)
print("Input mode:", input_provenance["mode"])
print("Geocoding acquired:", input_provenance["geocoding_retrieved_utc"])
data = {
"address": ["Centraal Station", "VU Hoofdgebouw", "Amsterdam Business School",
"University of Amsterdam", "Johan Cruijff ArenA", "Ziggo Dome"],
"color": ["blue", "black", "orange", "green", "red", "purple"],
}
df = pd.DataFrame(data)
df["city"], df["country"] = "Amsterdam", "NL"
geocoding_record = json.loads((inputs_dir / "geocoding.json").read_text(encoding="utf8"))
features = geocoding_record["results"]
if len(features) != len(df):
raise ValueError("The graph input bundle must include all six place responses.")
df["lon"] = [feature["geometry"]["coordinates"][0] for feature in features]
df["lat"] = [feature["geometry"]["coordinates"][1] for feature in features]
display(style_table(df).format({"lat": "{:.6f}", "lon": "{:.6f}"}))
(evidence / "geocoding.json").write_bytes((inputs_dir / "geocoding.json").read_bytes())
Map = folium.Map(location=(52.37, 4.89), zoom_start=13)
for _, row in df.iterrows():
folium.Marker((row.lat, row.lon), icon=folium.Icon(color=row.color),
tooltip=row.address).add_to(Map)
Map
Input mode: recorded
Geocoding acquired: 2026-09-24T19:43:29.223917+00:00
| address | color | city | country | lon | lat | |
|---|---|---|---|---|---|---|
| 0 | Centraal Station | blue | Amsterdam | NL | 4.900871 | 52.379045 |
| 1 | VU Hoofdgebouw | black | Amsterdam | NL | 4.866900 | 52.334180 |
| 2 | Amsterdam Business School | orange | Amsterdam | NL | 4.911718 | 52.365107 |
| 3 | University of Amsterdam | green | Amsterdam | NL | 4.889804 | 52.368133 |
| 4 | Johan Cruijff ArenA | red | Amsterdam | NL | 4.942843 | 52.314353 |
| 5 | Ziggo Dome | purple | Amsterdam | NL | 4.937709 | 52.313753 |
Geocoding returns coordinates for a description, not necessarily the entrance we had in mind. Inspect the map and the returned locations. The default reuses the six Photon responses from the recorded input bundle. The printed input mode and acquisition date identify whether this run used those responses or a successful live refresh. The session export preserves that provenance.
ox.settings.log_console = False
ox.settings.use_cache = True
ox.settings.cache_folder = evidence / ".cache"
ox.settings.requests_timeout = 30
# Both input modes supply a complete graph/geocoding bundle.
G_walk = ox.load_graphml(inputs_dir / "amsterdam-walk.graphml.gz")
(evidence / "amsterdam-walk.graphml.gz").write_bytes(
(inputs_dir / "amsterdam-walk.graphml.gz").read_bytes()
)
print(f"G_walk has {G_walk.number_of_nodes():,} nodes and {G_walk.number_of_edges():,} arcs.")
print("Graph SHA-256:", input_provenance["graph_sha256"])
G_walk has 51,795 nodes and 143,968 arcs.
Graph SHA-256: 21e5c176f32a90127110e8c9d0049271df5614620629cd5fb17f678a4eb3760b
The graph comes from OpenStreetMap contributors (ODbL). Its municipal boundary and retained connected component can affect coverage. A route is shortest within this graph; it is not a guarantee about temporary closures or walking access at a venue.
# Visualize the graph
node_positions = {
node: (data["x"], data["y"]) for node, data in G_walk.nodes(data=True)
}
nx.draw_networkx_nodes(G_walk, pos=node_positions, node_size=1)
plt.show()
def style_locations(table):
return style_table(table).format({"lat": "{:.6f}", "lon": "{:.6f}",
"node": "{:.0f}", "snap_m": "{:.3f}"})
df["node"], df["snap_m"] = ox.distance.nearest_nodes(
G_walk, df.lon, df.lat, return_dist=True
)
display(style_locations(df[["address", "lat", "lon", "node", "snap_m"]]))
shared = df[df.duplicated("node", keep=False)]
if not shared.empty:
print("Distinct places assigned to the same graph node:")
display(style_locations(shared[["address", "node", "snap_m"]]))
# A configurable screening threshold; a short snap alone does not prove coverage.
MAX_SNAP_DISTANCE_M = 1000.0
if not np.isfinite(MAX_SNAP_DISTANCE_M) or MAX_SNAP_DISTANCE_M < 0:
raise ValueError("MAX_SNAP_DISTANCE_M must be finite and nonnegative.")
df["distant_snap"] = df["snap_m"] > MAX_SNAP_DISTANCE_M
print(f"Southernmost graph node: {min(a['y'] for _, a in G_walk.nodes(data=True)):.5f} N")
if df["distant_snap"].any():
print(f"Coverage warning: snaps exceed {MAX_SNAP_DISTANCE_M:.0f} m. "
"Routes involving these places are unsuitable as walking-distance estimates.")
display(style_locations(df.loc[df.distant_snap, ["address", "snap_m"]]))
| address | lat | lon | node | snap_m | |
|---|---|---|---|---|---|
| 0 | Centraal Station | 52.379045 | 4.900871 | 5629071965 | 11.179 |
| 1 | VU Hoofdgebouw | 52.334180 | 4.866900 | 294648098 | 61.952 |
| 2 | Amsterdam Business School | 52.365107 | 4.911718 | 46356661 | 43.880 |
| 3 | University of Amsterdam | 52.368133 | 4.889804 | 13285036240 | 55.953 |
| 4 | Johan Cruijff ArenA | 52.314353 | 4.942843 | 2680554148 | 2329.508 |
| 5 | Ziggo Dome | 52.313753 | 4.937709 | 2680554148 | 2016.547 |
Distinct places assigned to the same graph node:
| address | node | snap_m | |
|---|---|---|---|
| 4 | Johan Cruijff ArenA | 2680554148 | 2329.508 |
| 5 | Ziggo Dome | 2680554148 | 2016.547 |
Southernmost graph node: 52.31935 N
Coverage warning: snaps exceed 1000 m. Routes involving these places are unsuitable as walking-distance estimates.
| address | snap_m | |
|---|---|---|
| 4 | Johan Cruijff ArenA | 2329.508 |
| 5 | Ziggo Dome | 2016.547 |
Shortest paths and contraction hierarchies#
The routing problem is a weighted shortest-path problem. Let \(G=(V,A)\) be the walking network, with nonnegative arc lengths \(\ell_{uv}\). For snapped origin \(s\) and destination \(t\), let \(\mathcal P_{st}\) contain the paths from \(s\) to \(t\) that follow the allowed arcs and their directions. We seek
A feasible route stays on the supplied walking network and respects its allowed directions; if no such path exists, the destination is unreachable. Repeating this problem for each selected origin and destination gives the distance matrix. The snapping segments are outside this network objective; the access-distance extension below treats them separately. The timing experiment uses mutually reachable nodes and identical integer-millimetre arc costs for both implementations.
For a single route Dijkstra’s algorithm is straightforward. For many routes,
Pandana invests in contraction-hierarchy preprocessing. We record that cost
separately from the subsequent queries. The walking graph already contains its
allowed directions, so twoway=False avoids adding new reverse arcs.
construction_start = perf_counter()
nodes = pd.DataFrame.from_dict(dict(G.nodes(data=True)), orient="index")
edges = nx.to_pandas_edgelist(G)
network_start = perf_counter()
network = pandana.Network(nodes["x"], nodes["y"], edges["source"], edges["target"],
edges[["length"]], twoway=False)
network_construction_s = perf_counter() - network_start
preprocessing_s = perf_counter() - construction_start
print(f"DataFrame preparation + network construction: {preprocessing_s:.3f} s")
print(f"Network constructor (including hierarchy): {network_construction_s:.3f} s")
DataFrame preparation + network construction: 6.234 s
Network constructor (including hierarchy): 4.392 s
df["pandana"] = df["node"]
display(style_table(df[["address", "node", "pandana", "snap_m"]]))
| address | node | pandana | snap_m | |
|---|---|---|---|---|
| 0 | Centraal Station | 5629071965 | 5629071965 | 11.179360 |
| 1 | VU Hoofdgebouw | 294648098 | 294648098 | 61.952308 |
| 2 | Amsterdam Business School | 46356661 | 46356661 | 43.879517 |
| 3 | University of Amsterdam | 13285036240 | 13285036240 | 55.952679 |
| 4 | Johan Cruijff ArenA | 2680554148 | 2680554148 | 2329.507539 |
| 5 | Ziggo Dome | 2680554148 | 2680554148 | 2016.546853 |
origin, destination = df.loc[2, "node"], df.loc[3, "node"]
path_pandana = list(network.shortest_path(origin, destination))
path_nx = nx.shortest_path(G, origin, destination, weight="length")
Two optimal routes need not visit identical nodes. Check their endpoints, arcs, and costs rather than the set difference of their nodes. Pandana represents costs as integer thousandths of the supplied unit: here it truncates each arc to whole millimetres. For exact comparison we give NetworkX those same integer costs. This numerical convention does not imply millimetre accuracy of the map.
for _, _, attributes in G.edges(data=True):
attributes["length_mm"] = int(attributes["length"] * 1000)
assert path_pandana[0] == origin and path_pandana[-1] == destination
assert all(G.has_edge(u, v) for u, v in zip(path_pandana[:-1], path_pandana[1:]))
reference = nx.shortest_path_length(G, origin, destination, weight="length_mm") / 1000
assert np.isclose(nx.path_weight(G, path_pandana, "length_mm") / 1000,
reference, atol=1e-6, rtol=0)
print(f"NetworkX, original lengths: {nx.path_weight(G, path_nx, 'length'):.3f} m")
print(f"Pandana, integer-mm costs: {network.shortest_path_length(origin, destination):.3f} m")
NetworkX, original lengths: 1883.315 m
Pandana, integer-mm costs: 1883.298 m
# Reachability and geographic coverage are separate checks.
ids = df.node.to_numpy(dtype=np.int64)
if not all(nx.has_path(G, int(i), int(j)) for i in ids for j in ids):
raise ValueError("Some selected graph nodes are unreachable; inspect network components.")
distance_labels = df.address + df.distant_snap.map({True: " [coverage warning]", False: ""})
origs = np.repeat(ids, len(ids))
dests = np.tile(ids, len(ids))
distances = np.asarray(network.shortest_path_lengths(origs, dests)).reshape(len(ids), len(ids))
display(style_table(pd.DataFrame(distances, index=distance_labels, columns=distance_labels).rename_axis(index=None, columns=None).round(1)).format("{:.1f}").set_caption("Graph-node distances (m); flagged places lack nearby network coverage"))
| Centraal Station | VU Hoofdgebouw | Amsterdam Business School | University of Amsterdam | Johan Cruijff ArenA [coverage warning] | Ziggo Dome [coverage warning] | |
|---|---|---|---|---|---|---|
| Centraal Station | 0.0 | 6876.2 | 2288.7 | 1588.8 | 8689.6 | 8689.6 |
| VU Hoofdgebouw | 6876.2 | 0.0 | 6208.5 | 5425.1 | 7564.6 | 7564.6 |
| Amsterdam Business School | 2288.7 | 6208.5 | 0.0 | 1883.3 | 7376.6 | 7376.6 |
| University of Amsterdam | 1588.8 | 5425.1 | 1883.3 | 0.0 | 7543.5 | 7543.5 |
| Johan Cruijff ArenA [coverage warning] | 8689.6 | 7564.6 | 7376.6 | 7543.5 | 0.0 | 0.0 |
| Ziggo Dome [coverage warning] | 8689.6 | 7564.6 | 7376.6 | 7543.5 | 0.0 | 0.0 |
An access-distance extension#
One simple estimate adds the two straight-line snapping distances to the network distance: \(\widehat d_{ij}=s_i+d_{ij}^{\mathrm{network}}+s_j\) for distinct places. The distance from a place to itself remains zero. If two places share a node, the estimate includes travelling from each point to that common node.
These connectors are not verified walkable paths and can cross barriers. The estimate need not be the shortest real-world route, especially when the graph omits nearby streets. Comparing the two tables makes the effect of snapping visible rather than hiding it by replacing the example locations.
For the archived ArenA and Ziggo Dome coordinates, which are about 355 m apart, the two connectors sum to about 4.35 km. This demonstrates the coverage failure: it is not a usable walking-distance estimate between the venues. Rows and columns with distant snaps are explicitly flagged below.
access_estimate = distances + df.snap_m.to_numpy()[:, None] + df.snap_m.to_numpy()[None, :]
np.fill_diagonal(access_estimate, 0)
display(style_table(pd.DataFrame(access_estimate, index=distance_labels, columns=distance_labels).rename_axis(index=None, columns=None).round(1)).format("{:.1f}").set_caption("Connector illustration (m), not verified walking distances; flagged venue routes are unsuitable"))
| Centraal Station | VU Hoofdgebouw | Amsterdam Business School | University of Amsterdam | Johan Cruijff ArenA [coverage warning] | Ziggo Dome [coverage warning] | |
|---|---|---|---|---|---|---|
| Centraal Station | 0.0 | 6949.3 | 2343.7 | 1655.9 | 11030.3 | 10717.4 |
| VU Hoofdgebouw | 6949.3 | 0.0 | 6314.3 | 5543.0 | 9956.1 | 9643.1 |
| Amsterdam Business School | 2343.7 | 6314.3 | 0.0 | 1983.1 | 9750.0 | 9437.1 |
| University of Amsterdam | 1655.9 | 5543.0 | 1983.1 | 0.0 | 9928.9 | 9616.0 |
| Johan Cruijff ArenA [coverage warning] | 11030.3 | 9956.1 | 9750.0 | 9928.9 | 0.0 | 4346.1 |
| Ziggo Dome [coverage warning] | 10717.4 | 9643.1 | 9437.1 | 9616.0 | 4346.1 | 0.0 |
A batch of routes#
We first demonstrate a matrix for 50 sampled nodes. Sampling within one strongly connected component ensures that all selected pairs are reachable. Later we increase the matrix size systematically. Retrieve only a few full paths here; storing every path requires much more memory than storing the distance matrix.
component = np.array(sorted(max(nx.strongly_connected_components(G), key=len)), dtype=np.int64)
rng = np.random.default_rng(2025)
sample = rng.choice(component, min(50, len(component)), replace=False)
origs, dests = np.repeat(sample, len(sample)), np.tile(sample, len(sample))
table = pd.DataFrame(np.asarray(network.shortest_path_lengths(origs, dests)).reshape(len(sample), -1),
index=sample, columns=sample)
departure, arrival = np.unravel_index(np.argmax(table.to_numpy()), table.shape)
path_pandana = list(network.shortest_path(sample[departure], sample[arrival]))
print(f"Largest distance in this sample: {table.iloc[departure, arrival] / 1000:.3f} km")
paths = network.shortest_paths(origs[:5], dests[:5])
Largest distance in this sample: 24.141 km
street_names = []
for u, v in zip(paths[1][:-1], paths[1][1:]):
name = G[u][v].get("name", "Unnamed path")
if isinstance(name, list):
name = " / ".join(map(str, name))
if not street_names or name != street_names[-1]:
street_names.append(name)
print(" → ".join(street_names[:20]))
if len(street_names) > 20:
print(f"... {len(street_names) - 20} further named segments")
Unnamed path → Bornhout → Australiëhavenweg → Westpoortweg → Luvernes → Unnamed path → Sierenborch → Theemsweg → Seineweg → Naritaweg → Seineweg → Burgemeester De Vlugtlaan → Unnamed path → Burgemeester De Vlugtlaan → Bos en Lommerweg / Burgemeester De Vlugtlaan → Bos en Lommerweg → Bos en Lommerplein → Unnamed path → Hoofdweg → Unnamed path
... 31 further named segments
# These sampled endpoints already lie on the road network.
draw_route(paths[1])
How does the distance-matrix computation scale?#
The Amsterdam tables and timing plots below come from a fresh execution of the quick profile using the recorded inputs. They are new measurements, not a replay of the archived 24 September 2026 timing rows. The archive linked at the end preserves those earlier full-profile measurements separately.
The saved quick-profile measurements were recorded locally on 3 November 2026. Their timing rows, selected nodes and environment manifest are supplied separately from the September full run. The manifest identifies the same bundled graph by its checksum; it records the measurement session, while the graph and geocoding acquisition dates remain those of the recorded inputs.
The next cell records the machine and software for the current execution, including its CPU allocation. Colab hardware varies between sessions. The archived 24 September 2026 full-profile run was measured on an Intel Core Ultra 9 285HX, 24 cores / 24 logical processors, 63.46 GiB RAM, 64-bit Windows 11 build 26200, Python 3.14.6, Pandana 0.8, NetworkX 3.6.1, NumPy 2.5.3, pandas 3.0.5 and OSMnx 2.1.1. A new execution retains its own environment record; Windows measurements must not be relabelled as Colab. CPU affinity, power/frequency settings and background activity are not controlled.
The saved Amsterdam walking snapshot has 51,795 nodes and 142,338 directed arcs after keeping the shortest of parallel arcs (143,968 arcs before reduction). For a live download, use the graph counts printed by that execution. For each size \(n\), all methods compute the same \(n\times n\) matrix. NetworkX runs one single-source Dijkstra search per origin, exploring the reachable graph and then extracting the selected destinations. Pandana receives all \(n^2\) ordered pairs in one batch call. NetworkX’s standard Dijkstra routine is single-threaded, and our loop processes one origin at a time. We deliberately restrict Pandana’s construction and queries to one thread. The purpose is to showcase the different algorithmic implementations without parallelism as an additional source of speedup. The Python/C++ language difference remains.
Nested samples use seed 2025. The quick profile uses sizes 10, 25 and 50; the optional full profile uses 10, 25, 50, 100, 250, 500 and 1000. Each method has one warm-up and three measured repetitions per size. The order alternates across repetitions. Query wall time includes query-array preparation, distance computation, matrix conversion and temporary-object cleanup. Graph loading, hierarchy preprocessing, validation and plotting are excluded. Preprocessing is recorded separately, once. All methods use the same integer-millimetre arc costs; every matrix is checked outside the timer.
NetworkX uses one sequential Python worker; Pandana also uses compiled C++. Even the single-thread comparison therefore measures these implementations, not the algorithmic effect of contraction hierarchies in isolation. The sampling pool is fixed independently of the tested sizes; no matrix exceeds 1,000 locations.
hardware = {
"cpu": platform.processor() or "unavailable",
"physical_cores": psutil.cpu_count(logical=False),
"logical_processors": psutil.cpu_count(logical=True),
"ram_GiB": round(psutil.virtual_memory().total / 2**30, 2),
"os": platform.platform(), "architecture": platform.machine(),
}
if platform.system() == "Windows":
hardware["cpu"] = subprocess.check_output([
"powershell.exe", "-NoProfile", "-Command",
"(Get-CimInstance Win32_Processor).Name",
], text=True).strip()
if platform.system() == "Linux" and Path("/proc/cpuinfo").exists():
for line in Path("/proc/cpuinfo").read_text().splitlines():
if line.startswith("model name"):
hardware["cpu"] = line.split(":", 1)[1].strip()
break
try:
affinity_cpus = len(os.sched_getaffinity(0))
except AttributeError:
affinity_cpus = len(psutil.Process().cpu_affinity()) if hasattr(psutil.Process(), "cpu_affinity") else psutil.cpu_count()
cpu_quota = None
quota_file = Path("/sys/fs/cgroup/cpu.max")
if quota_file.exists():
quota, period = quota_file.read_text().split()
if quota != "max":
cpu_quota = int(quota) / int(period)
hardware["CPUs_available_by_affinity"] = affinity_cpus
hardware["cgroup_CPU_quota"] = cpu_quota
hardware["runtime_interface"] = "Google Colab" if IN_COLAB else "Local Python"
hardware["Colab_release_tag"] = os.environ.get("COLAB_RELEASE_TAG", "not reported")
environment = {
**hardware, "python": platform.python_version(),
"packages": {p: version(p) for p in ["pandana", "networkx", "numpy", "pandas", "osmnx"]},
"pandana_OpenMP_threads": PANDANA_THREADS, "NetworkX_workers": 1,
"seed": 2025, "repetitions": 3, "warmups_per_size": 1,
"benchmark_profile": "full" if RUN_FULL_BENCHMARK else "quick",
"input_provenance": input_provenance,
"network_nodes": len(G), "directed_arcs": G.number_of_edges(),
"graph_sha256": hashlib.sha256((evidence / "amsterdam-walk.graphml.gz").read_bytes()).hexdigest(),
"run_utc": datetime.now(timezone.utc).isoformat(),
"preprocessing_s": preprocessing_s, "network_constructor_s": network_construction_s,
}
display(style_table(pd.Series({**hardware, "Python": platform.python_version(),
"Pandana threads": PANDANA_THREADS,
**environment["packages"]}, name="experimental environment").to_frame()))
(evidence / "single-thread-environment.json").write_text(json.dumps(environment, indent=2), encoding="utf-8");
| experimental environment | |
|---|---|
| cpu | 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz |
| physical_cores | 4 |
| logical_processors | 8 |
| ram_GiB | 63.730000 |
| os | Windows-11-10.0.26300-SP0 |
| architecture | AMD64 |
| CPUs_available_by_affinity | 8 |
| cgroup_CPU_quota | None |
| runtime_interface | Local Python |
| Colab_release_tag | not reported |
| Python | 3.12.14 |
| Pandana threads | 1 |
| pandana | 0.8 |
| networkx | 3.6.1 |
| numpy | 2.5.3 |
| pandas | 3.0.5 |
| osmnx | 2.1.1 |
def dijkstra_matrix(selected):
matrix = np.empty((len(selected), len(selected)), dtype=float)
for i, origin in enumerate(selected):
lengths = nx.single_source_dijkstra_path_length(G, int(origin), weight="length_mm")
matrix[i] = [lengths[int(d)] / 1000 for d in selected]
del lengths
return matrix, {}
setup_times = {1: preprocessing_s}
def pandana_matrix(selected):
start = perf_counter()
origins = np.repeat(selected, len(selected))
destinations = np.tile(selected, len(selected))
prepared = perf_counter()
values = network.shortest_path_lengths(origins, destinations)
queried = perf_counter()
matrix = np.asarray(values).reshape(len(selected), -1)
converted = perf_counter()
del values, origins, destinations
cleaned = perf_counter()
return matrix, {"arrays_s": prepared-start, "batch_call_s": queried-prepared,
"conversion_s": converted-queried, "cleanup_s": cleaned-converted}
requested_sizes = [10, 25, 50, 100, 250, 500, 1000] if RUN_FULL_BENCHMARK else [10, 25, 50]
sizes = [n for n in requested_sizes if n <= len(component)]
ordered_sample = np.random.default_rng(2025).choice(
component, min(5000, len(component)), replace=False)[:max(sizes)]
methods = [("NetworkX (1 thread)", dijkstra_matrix), ("Pandana (1 thread)", pandana_matrix)]
np.save(evidence / "single-thread-selected.npy", ordered_sample)
environment["selected_sha256"] = hashlib.sha256(ordered_sample.tobytes()).hexdigest()
(evidence / "single-thread-environment.json").write_text(json.dumps(environment, indent=2), encoding="utf-8")
measurements = []
for n in sizes:
selected = ordered_sample[:n]
expected, _ = dijkstra_matrix(selected) # NetworkX warm-up and reference
warmup, _ = pandana_matrix(selected)
assert np.allclose(warmup, expected, atol=1e-6, rtol=0)
del warmup
for repetition in range(3):
order = methods if repetition % 2 == 0 else methods[::-1]
for method, function in order:
pools = [p for p in threadpool_info() if p["user_api"] == "openmp"]
assert pools and all(p["num_threads"] == 1 for p in pools)
gc.collect()
cpu_before = psutil.cpu_percent(interval=0.1)
start = perf_counter()
result, phases = function(selected)
seconds = perf_counter() - start
error = float(np.max(np.abs(result - expected)))
assert np.isfinite(result).all() and error <= 1e-6
measurements.append({"n": n, "method": method, "threads": 1,
"repetition": repetition+1, "seconds": seconds,
"system_cpu_before_pct": cpu_before,
"max_error_m": error, **phases})
del result
pd.DataFrame(measurements).to_csv(evidence / "single-thread-repetitions.csv", index=False)
print(f"Validated {n} x {n}: NetworkX and Pandana, both single-threaded.", flush=True)
del expected
timings = pd.DataFrame(measurements)
(evidence / "single-thread-setup.json").write_text(json.dumps(setup_times, indent=2));
Validated 10 x 10: NetworkX and Pandana, both single-threaded.
Validated 25 x 25: NetworkX and Pandana, both single-threaded.
Validated 50 x 50: NetworkX and Pandana, both single-threaded.
stats = timings.groupby(["n", "method"]).seconds.agg(["median", "min", "max"])
summary = stats["median"].unstack("method")
display(style_table(summary).format("{:.3f}").set_caption("Amsterdam: median query seconds, three repetitions"))
display(style_table(pd.Series(setup_times, name="Preprocessing seconds").rename_axis("Pandana threads").to_frame()))
summary.to_csv(evidence / "single-thread-summary.csv")
fig, axes = plt.subplots(1, 2, figsize=(10, 4), layout="constrained")
for method, color in [("NetworkX (1 thread)", "#b13b32"),
("Pandana (1 thread)", "#176a8a")]:
values = stats.xs(method, level="method")
pairs = values.index.to_numpy()**2
axes[0].plot(pairs, values["median"], "o-", label=method, color=color)
axes[0].fill_between(pairs, values["min"], values["max"], color=color, alpha=.18)
for threads, color in [(1, "#176a8a")]:
values = summary[f"Pandana ({threads} thread{'s' if threads > 1 else ''})"]
axes[1].plot(summary.index, summary["NetworkX (1 thread)"]/values,
"o-", color=color, label=f"{threads} thread(s): query only")
axes[1].plot(summary.index, summary["NetworkX (1 thread)"]/(values+setup_times[threads]),
"--", color=color, label=f"{threads} thread(s): with setup")
axes[0].set(xscale="log", yscale="log", xlabel="Origin-destination pairs (n²)",
ylabel="Query wall time (s)", title="Amsterdam; bands show min–max")
axes[1].set(xscale="log", yscale="log", xlabel="Locations (n)", ylabel="NetworkX / Pandana time",
title="Effect of charging preprocessing to one matrix")
axes[1].set_xticks(summary.index.to_numpy(), labels=[str(n) for n in summary.index])
axes[1].axhline(1, color="gray", linewidth=.8)
for ax in axes: ax.legend(fontsize=8)
plt.show()
| method | NetworkX (1 thread) | Pandana (1 thread) |
|---|---|---|
| n | ||
| 10 | 2.195 | 0.016 |
| 25 | 5.803 | 0.026 |
| 50 | 11.718 | 0.060 |
| Preprocessing seconds | |
|---|---|
| Pandana threads | |
| 1 | 6.234361 |
Interpreting the scaling experiment#
The query-only comparison asks how quickly each method answers a matrix request once its network is ready. The dashed speedup curve charges the complete measured Pandana preprocessing cost to one matrix; this cost is paid once when the same network serves repeated requests. The spread over repetitions indicates observed timing variability, not a confidence interval. Preprocessing is measured once.
The quick profile is a short illustration, not a measurement of the full scaling range. Enable the full profile to examine larger matrices on this machine. The earlier archived full run is separate historical evidence; its timings should not be interpreted as results of this execution. On a fixed graph, NetworkX performs one full search per origin while Pandana answers all \(n^2\) pairs. Both also incur quadratic matrix/output work. This is not a general linear-versus- quadratic complexity comparison.
Increasing \(n\) increases the number of pairs quadratically. Even fast routing does not remove the cost of representing the result: the final float64 matrix alone needs \(8n^2\) bytes. Pandana additionally constructs origin and destination arrays and returns a Python list before conversion, so the final matrix size is not peak memory usage. NetworkX temporarily holds a one-to-many distance dictionary for each row. Neither method materializes all full paths in this experiment.
Use the measured curves to assess when preprocessing pays off on this graph. Do not extrapolate the observed speedup to arbitrarily large matrices or attribute all of it to contraction hierarchies alone: language, implementation overhead, and the amount of graph explored also differ. Try increasing the sample sizes if memory and runtime permit. Keeping both methods single-threaded makes the routing-implementation comparison easier to interpret.
Where does Pandana spend the query time?#
The same measured calls above record three phases: preparing ordered-pair arrays, the public Pandana batch call, and converting the returned list into a NumPy matrix. The batch call includes the Python wrapper, validation and native routing; it is not a pure contraction-hierarchy kernel timer. Temporary-object cleanup is recorded separately, and the end-to-end query timer includes it. Preprocessing is outside all query timers. Phase medians need not sum to the median total.
The following table and plot use one thread throughout. They are a breakdown of the current measurements, rather than a second experiment with different timings.
phase_summary = timings[timings.method.str.startswith("Pandana")].groupby(["n", "threads"])[
["arrays_s", "batch_call_s", "conversion_s", "cleanup_s", "seconds"]].median()
display(style_table(phase_summary).format("{:.6f}").set_caption("Amsterdam: median phase seconds, three repetitions"))
phase_summary.to_csv(evidence / "single-thread-phases.csv")
| arrays_s | batch_call_s | conversion_s | cleanup_s | seconds | ||
|---|---|---|---|---|---|---|
| n | threads | |||||
| 10 | 1 | 0.000085 | 0.015599 | 0.000017 | 0.000002 | 0.015743 |
| 25 | 1 | 0.000090 | 0.025893 | 0.000062 | 0.000008 | 0.026076 |
| 50 | 1 | 0.000106 | 0.059836 | 0.000176 | 0.000027 | 0.060156 |
fig, ax = plt.subplots(figsize=(8, 4), layout="constrained")
parts = phase_summary.xs(1, level="threads")
for column, label in [("arrays_s", "Prepare pair arrays"),
("batch_call_s", "Pandana batch call"), ("conversion_s", "Convert to NumPy")]:
ax.plot(parts.index, parts[column], "o-", label=label)
ax.set(xscale="log", yscale="log", xlabel="Locations (n)", ylabel="Median phase time (s)",
title="Pandana query phases: one thread")
ax.legend()
plt.show()
A larger graph: The whole Netherlands#
The following recorded Windows results were obtained on 24 September 2026,
outside Colab, on an Intel Core Ultra 9 285HX (24 cores / 24 logical processors),
63.46 GiB RAM, Windows 11 build 26200, with Python 3.14.6, Pandana 0.8,
NetworkX 3.6.1, NumPy 2.5.3, pandas 3.0.5, and OSMnx 2.1.1. They use the
Geofabrik Netherlands PBF dated 23 September 2026.
That historical download URL returned HTTP 404 when checked on 27 October 2026.
Its SHA-256 is preserved in the archived source.json, but the national PBF,
derived graph arrays and reference matrices are not currently available in a
public archive. The national figures therefore remain recorded results, not
an exactly reproducible public-input experiment. Substituting today’s extract
would create a new experiment. The separate Amsterdam graph and inputs are
available in the saved archive linked at the end of this notebook.
The source is OpenStreetMap contributors, ODbL 1.0. The walking graph has
3,333,231 nodes, 8,454,986 directed arcs, and 33,821 connected components.
We retain the entire graph; seed 2025 selects 1,000 distinct endpoints in the
largest component of 3,143,652 nodes (94.31%). All one million selected pairs
are reachable. These are 1,000-by-1,000 matrices on the national graph, not an
all-nodes distance matrix.
The walking filter matches OSMnx’s walking filter; arcs are bidirectional, and ferries and relation turn restrictions are excluded. The streaming PBF builder retains way endpoints and shared junctions, sums lengths along intermediate geometry, and keeps the shortest parallel edge. Unlike Amsterdam’s OSMnx simplification, it does not remove degree-two nodes across way boundaries. Consequently this is not a controlled graph-size-only comparison.
For this larger graph as well, we deliberately restrict Pandana to one thread to match the single-threaded NetworkX Dijkstra routine and sequential origin loop, keeping the focus on the routing implementations rather than parallelism. NetworkX completed one cold, single-threaded run, with no warm-up or restart. Its query time sums the 1,000 row timings, including Dijkstra, destination extraction, row-array creation and dictionary cleanup; checkpoint writes, validation and final matrix assembly are excluded. Pandana was measured in three sequential, independent repetitions. Each repetition started a new Python process, built its own hierarchy with one OpenMP thread, and performed one warm-up followed by one measured matrix. Each repetition finished before the next started; the reported median is not a concurrent multi-process result. Its query timer includes pair arrays, the public batch call, matrix conversion and cleanup. Graph loading and preprocessing are excluded for both methods. NetworkX setup is graph construction; Pandana setup additionally builds the contraction hierarchy. Pandana setup and query entries are separate medians. The library runs were not interleaved, and background load was not controlled.
Every result matched an independent SciPy 1.18.1 integer-millimetre Dijkstra reference to within \(10^{-6}\) m; the observed maximum difference was zero. The reference was generated with four validation workers outside the benchmark timers. The cells below display these saved measurements; rerunning the notebook does not download the national PBF or repeat the two-hour NetworkX experiment.
The archived experiment code and reproduction notes identify the exact scripts and original result metadata from this Windows run. The Pandana launcher executes the repetitions sequentially; its worker defines the timing boundaries. The NetworkX worker implements the sequential row-by-row measurement and checkpointing. The Amsterdam experiment code is included directly in this notebook.
national = pd.DataFrame([{'method': 'NetworkX (1 thread)', 'query_s': 6572.529253199969, 'min_s': 6572.529253199969, 'max_s': 6572.529253199969, 'setup_s': 14.882086199999321, 'runs': 1}, {'method': 'Pandana (1 thread)', 'query_s': 17.8049972000008, 'min_s': 17.766781000000265, 'max_s': 18.27456390000043, 'setup_s': 114.2332279000002, 'runs': 3}]).set_index("method")
display(style_table(national).format({"query_s": "{:.3f}", "min_s": "{:.3f}", "max_s": "{:.3f}", "setup_s": "{:.3f}", "runs": "{:.0f}"}).set_caption("Netherlands: recorded seconds; NetworkX one run, Pandana median of three"))
fig, ax = plt.subplots(figsize=(8, 3.8), layout="constrained")
ax.barh(national.index, national.query_s, color=["#b13b32", "#176a8a"])
ax.set(xscale="log", xlabel="Query wall time (s; logarithmic scale)", title="Netherlands: 1,000 × 1,000 distances")
for i, value in enumerate(national.query_s):
ax.text(value*1.1, i, f"{value:,.2f} s", va="center")
ax.set_xlim(right=national.query_s.max()*6)
plt.show()
| query_s | min_s | max_s | setup_s | runs | |
|---|---|---|---|---|---|
| method | |||||
| NetworkX (1 thread) | 6572.529 | 6572.529 | 6572.529 | 14.882 | 1 |
| Pandana (1 thread) | 17.805 | 17.767 | 18.275 | 114.233 | 3 |
With one thread for both implementations, NetworkX took 6,572.53 seconds (109.54 minutes) and Pandana took 17.80 seconds, after preprocessing. The query-time ratio is approximately 369×; it describes these runs, not a confidence interval. Pandana’s single-thread preprocessing median was 114.23 seconds. This comparison showcases the routing implementations; Python versus compiled C++ remains an additional source of the difference.
Save this session’s results#
This cell packages the current Amsterdam timing records, selected nodes, graph checksum, graph snapshot, geocoding responses, software and hardware description, plus the already recorded national table with its separate Windows provenance. It does not claim a new national run. In Colab it downloads one ZIP file. Download the executed notebook from Colab’s File menu as well, so its plots and outputs are preserved with the measurements.
The saved local inputs and results form the archive of the 24 September 2026 full-profile Windows run of the Amsterdam example, together with the recorded national table. They do not include the national graph or its reference matrices. New executions generate their own archive below, recording both the input mode and the quick/full benchmark profile. The source data are © OpenStreetMap contributors, available under ODbL 1.0. The recorded graph’s known venue-coverage limitation is retained and flagged; offline availability does not make those venue routes suitable walking estimates.
assert environment["run_utc"] and len(timings) == 6 * len(sizes)
assert timings.max_error_m.max() <= 1e-6
session_stamp = environment["run_utc"].replace(":", "").replace("+", "_")
session_folder = evidence / ("session-" + session_stamp)
session_folder.mkdir(exist_ok=True)
timings.to_csv(session_folder / "amsterdam-repetitions.csv", index=False)
summary.to_csv(session_folder / "amsterdam-summary.csv")
phase_summary.to_csv(session_folder / "amsterdam-phases.csv")
np.save(session_folder / "amsterdam-selected.npy", ordered_sample)
national.to_csv(session_folder / "netherlands-recorded-windows.csv")
manifest = {
"amsterdam_environment": environment,
"input_provenance": input_provenance,
"data_attribution": "OpenStreetMap contributors, ODbL 1.0; https://www.openstreetmap.org/copyright",
"amsterdam_preprocessing_s": setup_times,
"exported_utc": datetime.now(timezone.utc).isoformat(),
"validation": {"measured_matrices": len(timings), "max_error_m": float(timings.max_error_m.max())},
"netherlands_provenance": {"execution": "Previously recorded standalone Windows runs, not this session",
"date": "2026-09-24", "cpu": "Intel Core Ultra 9 285HX", "physical_cores": 24,
"logical_processors": 24, "ram_GiB": 63.46, "os": "Windows 11 build 26200",
"python": "3.14.6", "pandana": "0.8", "networkx": "3.6.1",
"numpy": "2.5.3", "pandas": "3.0.5", "osmnx": "2.1.1"},
}
(session_folder / "manifest.json").write_text(json.dumps(manifest, indent=2), encoding="utf-8")
archive = evidence / (session_folder.name + ".zip")
with zipfile.ZipFile(archive, "w", compression=zipfile.ZIP_DEFLATED) as bundle:
for file in sorted(session_folder.iterdir()):
bundle.write(file, file.name)
for name in ["amsterdam-walk.graphml.gz", "geocoding.json"]:
bundle.write(evidence / name, name)
print(f"Saved {len(timings)} validated Amsterdam measurements and their provenance to {archive}")
if IN_COLAB:
files.download(str(archive))
else:
print("Open this ZIP from your notebook file browser; the published input archive is linked above.")
Saved 18 validated Amsterdam measurements and their provenance to data\shortest-path-live\session-2026-11-03T120547.540080_0000.zip
Open this ZIP from your notebook file browser; the published input archive is linked above.
Why graph size matters most#
The most important scalability benefit is handling much larger road graphs, especially for long-distance routes. Dijkstra repeatedly labels nodes across the reachable graph. Contraction hierarchies add distance-preserving shortcuts during preprocessing; at query time, forward and backward searches label nodes along upward edges toward higher-ranked vertices. Long routes can therefore traverse a much smaller high-level network, skipping most local detail. This progressively smaller search structure is the sense in which the graph becomes “sparser”; adding shortcuts need not reduce its stored edge count. See Geisberger’s description of contraction-hierarchy queries.
In the recorded Windows experiments, keeping the workload at one million pairs, the national single-thread query-time ratio is about 369×. NetworkX’s standard Dijkstra routine is single-threaded; we deliberately limit Pandana to one thread to showcase the algorithmic implementations without parallel speedup. Language, graph simplification and execution conditions still differ, so this is not a pure algorithm-only experiment; Colab results must be interpreted on their own machine.
From road distances to a traveling salesman tour#
The traveling salesman problem notebook shows how to use a matrix of pairwise distances to find a shortest tour visiting all selected locations and returning to the start. Combining the two notebooks lets us solve a TSP using real road-network distances: choose the locations, compute their distance matrix here, and supply it as the cost matrix to a TSP formulation there instead of using straight-line distances. Keep the same location ordering for the matrix rows and columns, and choose locations that are mutually reachable, with adequate coverage and reasonable snaps at every place. Exclude the flagged venue estimates above until their missing coverage is repaired. A directed road network can produce asymmetric distances, so use a directed TSP formulation unless the matrix is symmetric.