Routing Algorithms in Python: Dijkstra, A*, and Contraction Hierarchies

Picking the wrong shortest-path algorithm is the quietest way to blow a latency budget. A Dijkstra pass that expands uniformly across a metropolitan graph touches hundreds of thousands of nodes to answer a query that A* would settle in a few thousand; a contraction-hierarchy build that runs for twenty minutes is wasted preprocessing if your topology changes every hour. The three algorithms in this guide are not interchangeable — they trade optimality guarantees, preprocessing cost, and per-query node expansion against each other, and the correct choice is a function of graph size, query volume, and how often the edges move. This is the decision framework, grounded in runnable code: a hand-written priority-queue Dijkstra and A* that pull edges from an async Neo4j session, the Cypher and GDS equivalents, and the profiling discipline to know which one your workload actually needs. It sits under Network Routing Algorithms in Python.

Prerequisites

The examples run against a Neo4j 5.x instance with native point geometry and, for the library-backed path, the Graph Data Science plugin installed. The hand-written traversals need only the core async driver; the GDS variant needs the plugin and the graphdatascience client (or plain Cypher CALL gds.*).

pip install "neo4j>=5.18" "pytest>=8.0" "pytest-asyncio>=0.23"

The graph these algorithms traverse must already carry coordinates as native point values and keep a routing cost distinct from raw length — the modelling groundwork lives in node and edge spatial mapping, and the index that makes the endpoint lookups seek rather than scan comes from your spatial indexing strategy.

Core Concept & Mechanism

Every shortest-path algorithm in this family is a variation on the same loop: maintain a frontier of reachable nodes ordered by some key, pop the most promising one, relax its outgoing edges, and repeat until the destination is settled. What separates the algorithms is the ordering key and whether the graph is preprocessed first.

Dijkstra orders the frontier by $g(n)$, the confirmed lowest cost from the source to node $n$. Because it always settles the cheapest unsettled node next, it never has to revisit a settled node, and the first time it pops the destination that cost is provably optimal. The price is directional blindness: with nothing but $g$ to go on, the frontier grows as a roughly circular wavefront, so a query between two points 40 km apart on a continental graph will happily explore nodes 40 km in the wrong direction before the wave reaches the goal.

A* orders the frontier by $f(n) = g(n) + h(n)$, where $h(n)$ is a heuristic estimate of the remaining cost from $n$ to the destination. On a geographic graph the straight-line great-circle distance to the goal is a free, admissible heuristic — it never overestimates, because no road is shorter than the crow-flies path. Admissibility is exactly the property that preserves Dijkstra’s optimality guarantee while collapsing the frontier from a disk into a narrow ellipse aimed at the goal. The full derivation, and the scaling factor you need when edge weights are travel time rather than distance, is the subject of Implementing A* with a Haversine Heuristic in Python.

Contraction hierarchies (CH) attack the problem from the other side. Instead of a smarter per-query search, they precompute a node ordering by importance and insert shortcut edges that preserve shortest-path distances while letting a query skip whole chains of low-importance nodes. A bidirectional search over the contracted graph then answers point-to-point queries on a country-scale network in microseconds — but the preprocessing must be rebuilt when edge weights change, so CH fits static or slowly-changing topologies, not graphs under constant live edits.

$$\underbrace{g(n)}{\text{cost so far}} ;+; \underbrace{h(n)}{\substack{0 \text{ for Dijkstra} \ \text{admissible estimate for A*}}} ;=; f(n) \quad\text{(the frontier ordering key)}$$

Setting $h(n) = 0$ turns A* back into Dijkstra — they are the same algorithm with a different priority key, which is why a single priority-queue implementation can serve both.

Schema & Data Model

The traversal reads two things per edge: the neighbour it points to and the cost of crossing it. Store the cost as a precomputed scalar on the relationship so the search never recomputes geometry mid-loop, and keep it distinct from raw length so a time-based route and a distance-based route can share one graph.

CREATE CONSTRAINT junction_id_unique IF NOT EXISTS
FOR (j:Junction) REQUIRE j.id IS UNIQUE;

CREATE POINT INDEX junction_location IF NOT EXISTS
FOR (j:Junction) ON (j.location);

// Representative shape of the routable graph
// (:Junction {id, location: point({srid:4326, latitude, longitude})})
//   -[:ROAD {cost_s, length_m}]->
// (:Junction)

Here cost_s is the traversal weight the algorithms minimise (travel time in seconds), and length_m is retained separately for the heuristic and for reporting. Directionality is load-bearing: a one-way street is a single :ROAD relationship, a two-way street is two. The hand-written search below reads only outgoing relationships, so encoding direction in the edges is what keeps illegal manoeuvres out of the result.

Step-by-Step Implementation

We build one priority-queue search that runs as Dijkstra or A* depending on the heuristic passed to it, streaming a node’s neighbours from Neo4j on demand rather than loading the whole graph into memory. This is the shape you reach for when the graph is too large to materialise client-side but each query only touches a corridor of it.

1. Load a node’s outgoing edges asynchronously

The traversal calls this once per settled node. It returns each neighbour with the edge cost and the neighbour’s coordinates so an A* heuristic can be computed without a second round trip.

import asyncio
import heapq
import itertools
import math
from typing import Callable, Optional

from neo4j import AsyncGraphDatabase, AsyncSession

NEIGHBOUR_QUERY = """
MATCH (u:Junction {id: $node_id})-[r:ROAD]->(v:Junction)
RETURN v.id AS to_id,
       r.cost_s AS cost,
       v.location.latitude  AS lat,
       v.location.longitude AS lon
"""


async def expand(session: AsyncSession, node_id: str):
    """Yield (neighbour_id, edge_cost, lat, lon) for one node's out-edges."""
    result = await session.run(NEIGHBOUR_QUERY, node_id=node_id)
    return [
        (rec["to_id"], float(rec["cost"]), rec["lat"], rec["lon"])
        async for rec in result
    ]

2. Run the shared search loop

The heuristic is a plain callable. Pass lambda *_: 0.0 and the loop is Dijkstra; pass a great-circle estimate and it is A*. A monotonically increasing counter breaks ties on equal f, which keeps the heap from ever comparing node ids and makes the ordering deterministic.

async def shortest_path(
    session: AsyncSession,
    source: str,
    target: str,
    target_lat: float,
    target_lon: float,
    heuristic: Callable[[float, float], float],
) -> Optional[tuple[float, list[str]]]:
    """Uniform Dijkstra/A* search. heuristic(lat, lon) estimates cost to target."""
    counter = itertools.count()
    frontier: list[tuple[float, int, str]] = [(0.0, next(counter), source)]
    best_g: dict[str, float] = {source: 0.0}
    came_from: dict[str, str] = {}
    settled: set[str] = set()

    while frontier:
        f_score, _, node = heapq.heappop(frontier)
        if node in settled:
            continue
        if node == target:
            return best_g[node], _reconstruct(came_from, source, target)
        settled.add(node)

        for to_id, edge_cost, lat, lon in await expand(session, node):
            tentative = best_g[node] + edge_cost
            if tentative < best_g.get(to_id, math.inf):
                best_g[to_id] = tentative
                came_from[to_id] = node
                f = tentative + heuristic(lat, lon)
                heapq.heappush(frontier, (f, next(counter), to_id))
    return None


def _reconstruct(came_from: dict[str, str], source: str, target: str) -> list[str]:
    path, node = [target], target
    while node != source:
        node = came_from[node]
        path.append(node)
    path.reverse()
    return path

3. Drive it as either algorithm

The great-circle heuristic here estimates cost in the same unit as cost_s. If edges are travel time, the raw distance must be divided by the fastest plausible speed so the estimate never exceeds real travel time — that admissibility scaling is the crux of the dedicated A* page.

def great_circle_seconds(lat: float, lon: float, t_lat: float, t_lon: float,
                         max_mps: float = 33.3) -> float:
    """Admissible time estimate: crow-flies metres / top speed (m/s)."""
    R = 6_371_000.0
    p1, p2 = math.radians(lat), math.radians(t_lat)
    dphi = math.radians(t_lat - lat)
    dlam = math.radians(t_lon - lon)
    a = math.sin(dphi / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(dlam / 2) ** 2
    metres = 2 * R * math.asin(math.sqrt(a))
    return metres / max_mps


async def main():
    driver = AsyncGraphDatabase.driver(
        "neo4j://localhost:7687", auth=("neo4j", "secure-password"),
        max_connection_pool_size=20, connection_acquisition_timeout=15.0,
    )
    src, dst = "J-chi-0417", "J-chi-9920"
    dst_lat, dst_lon = 41.8500, -87.6500  # Chicago Loop
    try:
        async with driver.session(database="neo4j") as session:
            dijkstra = await shortest_path(
                session, src, dst, dst_lat, dst_lon, heuristic=lambda *_: 0.0)
            astar = await shortest_path(
                session, src, dst, dst_lat, dst_lon,
                heuristic=lambda lat, lon: great_circle_seconds(lat, lon, dst_lat, dst_lon))
        print(f"Dijkstra cost: {dijkstra[0]:.1f}s over {len(dijkstra[1])} nodes")
        print(f"A* cost:       {astar[0]:.1f}s over {len(astar[1])} nodes")
        assert abs(dijkstra[0] - astar[0]) < 1e-6  # same optimum, fewer expansions
    finally:
        await driver.close()


if __name__ == "__main__":
    asyncio.run(main())

Both calls return the identical optimal cost; A* simply pops fewer nodes off the frontier to get there. That equality is the assertion worth keeping in tests — if it ever fails, the heuristic has become inadmissible.

Query Patterns & Variants

The hand-written loop is one of three ways to run a shortest path against Neo4j. Match the tool to the workload.

Variant A — server-side shortestPath (unweighted). When you only need hop count and every edge is equal, the built-in is unbeatable because it never leaves the database. It minimises hops, not cost, so it is wrong the moment cost_s varies.

MATCH (s:Junction {id: $src}), (d:Junction {id: $dst})
MATCH p = shortestPath((s)-[:ROAD*..60]->(d))
RETURN [n IN nodes(p) | n.id] AS route, length(p) AS hops
// Hop-optimal only. Ignores cost_s entirely — never use it for weighted routing.

Variant B — GDS weighted Dijkstra. For weighted single-source or one-to-many routing at scale, project a subgraph once and let the library’s parallel Dijkstra run in-database over it. This is the production default when the corridor is large or you need cost-to-many-targets; the full projection lifecycle is covered in Weighted Dijkstra Routing with Neo4j GDS.

CALL gds.shortestPath.dijkstra.stream('road_proj', {
  sourceNode: gds.util.asNode($src_internal_id),
  targetNode: gds.util.asNode($dst_internal_id),
  relationshipWeightProperty: 'cost_s'
})
YIELD totalCost, nodeIds
RETURN totalCost, [id IN nodeIds | gds.util.asNode(id).id] AS route

Variant C — client-side A* (the loop above). When you want application-level control over the cost function — dynamic edge penalties, per-request road closures, custom tie-breaking — the hand-written search wins, because the logic lives in Python where you can change it per call without reprojecting a graph. It pays a network round trip per expanded node, so it is best on corridors bounded first by a distance filter rather than on unbounded continental searches.

The head-to-head numbers behind “when does GDS beat hand-written Cypher” live in Neo4j GDS vs Custom Cypher Routing.

Performance Tuning

The single number that predicts routing latency is nodes expanded, and PROFILE plus a settled-node counter are how you measure it. For the hand-written search, log len(settled) per query and watch the ratio between Dijkstra and A* on real origin-destination pairs — on a well-formed geographic graph A* should expand three to ten times fewer nodes; if the ratio is near one, the heuristic is not biasing the search and is almost certainly returning zero or a near-zero estimate.

• Confirm the endpoint lookup seeks. Both hand-written and GDS paths anchor on MATCH (:Junction {id}). Run PROFILE and confirm a NodeUniqueIndexSeek, never a NodeByLabelScan — a missing constraint on Junction.id turns every neighbour fetch into a full scan and dominates the whole query.
• Bound the corridor before you search. An unbounded client-side A* on a continental graph still round-trips per node. Pre-clip candidates with the box-then-distance predicate from distance filter query patterns so the frontier can never wander off the map.
• Batch expansion where possible. One round trip per settled node is the hand-written loop’s ceiling. If latency matters more than per-request cost-function flexibility, push the whole search into GDS where expansion never crosses the wire.
• Keep weights on the relationship. Recomputing cost_s from geometry inside the loop repeats trig on every edge of every candidate. Precompute it at ingestion and read it as a scalar.

Edge Cases & Gotchas

• Negative or zero edge weights. Dijkstra and A* both assume non-negative weights; a single negative cost_s (a mis-signed turn bonus, say) breaks the settled-node invariant and can return a wrong path with no error. Validate weight sign at ingestion, not at query time.
• Inadmissible heuristic silently loses optimality. If h(n) ever exceeds the true remaining cost, A* can settle the target on a sub-optimal path. The usual cause is a unit mismatch — a distance-metre heuristic against time-second edge weights. Scale the heuristic into the edge unit; the derivation is in the A* Haversine page.
• Disconnected components. A query between two nodes in different components expands the entire source component before returning None. Cap the search with a settled-node ceiling or a max-cost cutoff so an unreachable target fails fast.
• Stale contraction hierarchies. A CH built before an edge-weight update returns distances for the old graph. Treat the preprocessed structure as a cache keyed on a topology version and invalidate it on every write that touches routing weights.
• Directionality drift. Reading only out-edges is correct only if two-way streets are modelled as two relationships. A one-way edge stored bidirectionally lets the search cheat through it the wrong way.

Verification & Testing

Two properties are worth pinning in CI: A* returns the same cost as Dijkstra (optimality preserved), and both return the known-correct cost on a tiny hand-built graph. The second catches sign and unit bugs; the first catches an inadmissible heuristic sneaking in during a refactor.

import pytest
from neo4j import AsyncGraphDatabase

SEED = """
CREATE (a:Junction {id: 'A', location: point({srid:4326, latitude: 41.88, longitude: -87.64})})
CREATE (b:Junction {id: 'B', location: point({srid:4326, latitude: 41.89, longitude: -87.62})})
CREATE (c:Junction {id: 'C', location: point({srid:4326, latitude: 41.90, longitude: -87.63})})
CREATE (a)-[:ROAD {cost_s: 60.0, length_m: 1800.0}]->(b)
CREATE (b)-[:ROAD {cost_s: 40.0, length_m: 1200.0}]->(c)
CREATE (a)-[:ROAD {cost_s: 150.0, length_m: 4500.0}]->(c)
"""


@pytest.mark.asyncio
async def test_astar_matches_dijkstra_and_optimum():
    driver = AsyncGraphDatabase.driver("bolt://localhost:7687", auth=("neo4j", "test"))
    async with driver.session(database="neo4j") as s:
        await s.run("MATCH (n) DETACH DELETE n")
        await s.run(SEED)
        await s.run("CREATE CONSTRAINT jid IF NOT EXISTS "
                    "FOR (j:Junction) REQUIRE j.id IS UNIQUE")

        d = await shortest_path(s, "A", "C", 41.90, -87.63, heuristic=lambda *_: 0.0)
        a = await shortest_path(
            s, "A", "C", 41.90, -87.63,
            heuristic=lambda lat, lon: great_circle_seconds(lat, lon, 41.90, -87.63))

    assert d[0] == pytest.approx(100.0)      # A->B->C (60+40) beats direct 150
    assert d[1] == ["A", "B", "C"]
    assert a[0] == pytest.approx(d[0])       # A* preserves the optimum
    await driver.close()

Run the same assertion over a sample of real origin-destination pairs from production traffic, not just the toy graph — a heuristic that is admissible in one region can drift inadmissible where your speed assumptions break down (mountain passes, ferries, congestion-priced segments).

FAQ

When is Dijkstra actually the right choice over A*?

Use Dijkstra when you have no admissible heuristic, when edge weights have no geometric interpretation, or when you need cost from one source to many targets at once. A star only helps a single point-to-point query where a lower-bound estimate of remaining cost exists. For a one-to-many cost surface — every reachable node’s distance from a depot, for example — the heuristic has no single goal to bias toward, so plain Dijkstra is both correct and simpler.

Can I just set the A* heuristic to zero to reuse the same code?

Yes, and that is exactly why the implementation takes the heuristic as a parameter. Setting h(n) to zero makes f(n) equal to g(n), which is precisely Dijkstra. One priority-queue loop serves both algorithms; the only difference is the callable you pass in. This also makes testing easy, because you can assert that the zero-heuristic run and the geographic-heuristic run return the same cost.

Why does my A* return a slightly different route than Dijkstra?

The costs should be identical; if the node sequence differs it is almost always tie-breaking between two equal-cost optimal paths, which is harmless. If the total cost differs, the heuristic is inadmissible — it overestimated remaining cost somewhere, usually a unit mismatch between a distance heuristic and time-based edge weights. Scale the heuristic into the same unit as the edge weight and the costs will converge.

Should I load the whole graph into Python or query edges on demand?

On demand for large graphs, in memory for small stable ones. The streaming approach shown here keeps the client footprint flat and only touches the corridor a single query needs, at the cost of one round trip per expanded node. If the graph fits in memory and you run thousands of queries against it, load it once into an adjacency structure and skip the per-node round trips, or push the search into GDS where expansion never leaves the database.

How large does the graph need to be before contraction hierarchies pay off?

Preprocessing pays off when query volume on a stable graph is high enough to amortise the build. For occasional queries or a graph whose weights change frequently, the build cost never returns. As a rule of thumb, reach for contraction hierarchies on country- or continent-scale static road networks serving sustained high query rates, and stay with A star for interactive point-to-point routing on graphs that change often.

Related

• Implementing A* with a Haversine Heuristic in Python — the admissible heuristic and the time-versus-distance scaling factor in full.
• Weighted Dijkstra Routing with Neo4j GDS — projecting a subgraph and running library Dijkstra for one-to-many routing.
• Neo4j GDS vs Custom Cypher Routing — the benchmark that decides when to leave the database and when to stay.
• Distance Filter Query Patterns — bounding the corridor before a search so the frontier cannot wander.
• Spatial Graph Database Fundamentals for Python — the schema, indexing, and driver foundations these algorithms assume.

This guide is part of Network Routing Algorithms in Python.

For authoritative algorithm references, consult the Neo4j Graph Data Science path-finding documentation and the original contraction-hierarchies work by Geisberger et al..