07 - Algoritmi shortest path: Dijkstra, Bellman-Ford, Floyd-Warshall

Nell'articolo precedente la ricerca in ampiezza ha trovato il percorso con meno archi e un controesempio ha mostrato che non è il più economico: un collegamento diretto da 10 batte, per numero di archi, un giro da tre tappe da 1 ciascuna. Quando gli archi hanno un costo (chilometri, minuti, euro di pedaggio) la domanda diventa «qual è il percorso di costo minimo?», e la risposta dipende da una sola cosa: ci sono pesi negativi oppure no. Dijkstra è veloce ma non li ammette; Bellman-Ford li ammette e riconosce i cicli negativi; Floyd-Warshall risponde per tutte le coppie in una volta. Qui li implementi tutti e tre, li confronti con networkx su migliaia di grafi, e costruisci il caso in cui Dijkstra sbaglia. Ti servono le liste di adiacenza e la BFS dell'articolo 6.

Che cosa si calcola, e il gesto che li accomuna

Dato un grafo orientato con un peso su ogni arco, il costo di un percorso è la somma dei pesi dei suoi archi, e la distanza fra due vertici è il minimo di questi costi (infinita se non c'è percorso). Il gesto comune ai tre algoritmi è il rilassamento di un arco (u, v) di peso w: se arrivare a v passando da u costa meno della migliore stima che hai, cioè dist[u] + w < dist[v], aggiorni dist[v]. Le stime partono da infinito (0 sulla sorgente) e scendono soltanto; gli algoritmi differiscono per l'ordine in cui rilassano gli archi, ed è quell'ordine a decidere quando il risultato è garantito.

Il modulo sp.py, che apre il codice di tutto l'articolo, contiene le tre funzioni. Le leggerai a pezzi nelle sezioni che seguono; per ora basta sapere che ciascuna conta le proprie operazioni (inserimenti nello heap, passate sugli archi, passi del ciclo interno) perché la complessità si dichiara con un conteggio, non con i secondi.

Dijkstra: sempre l'estremo più vicino, finché i pesi sono non negativi

Dijkstra tiene uno heap di coppie (distanza, vertice) ed estrae sempre il vertice con la stima più piccola: a quel punto la sua stima è definitiva, perché ogni altro percorso passerebbe per vertici già più lontani e, con pesi non negativi, non potrebbe che aggiungere costo. Poi rilassa gli archi uscenti da quel vertice. Il paragone è un incendio che si propaga a velocità uguale al peso degli archi: ogni vertice si accende nell'ordine della sua distanza vera, e quello che si accende non si riaccende più.

In Python lo heap è heapq, e due dettagli del modulo contano. Non c'è un'operazione «abbassa la priorità»: si inserisce una voce nuova e le vecchie restano, quindi all'estrazione bisogna riconoscerle e scartarle, come suggerisce la documentazione di heapq con la tecnica della voce marcata. Inoltre le coppie si confrontano elemento per elemento, quindi a parità di distanza decide il vertice: con vertici numerici va bene, con oggetti non confrontabili servirebbe un contatore, come nell'esempio della stessa pagina.

La funzione dijkstra di sp.py ha due modalità. Con definitivo=True (la versione classica) un vertice estratto una volta viene messo in chiusi e non viene più espanso; con False si scartano solo le voci vecchie (d > dist[u]) e un vertice migliorato dopo l'estrazione viene espanso di nuovo. La differenza non conta con pesi non negativi e conta moltissimo con pesi negativi, come vedrai fra poco. Prima la prova di correttezza: 3.000 grafi orientati casuali con pesi interi da 0 a 20, confrontati con networkx (3.7, in un ambiente virtuale di prova), per entrambe le modalità. Poi i contatori su tre grafi grandi, e infine il caso peggiore.

import heapq

INF = float("inf")


def dijkstra(adj, s, definitivo=True):
    """adj[u] = lista di (v, peso). Ritorna distanze e contatori di operazioni.

    definitivo=True: un vertice estratto dallo heap non viene mai riaperto
    (la versione classica, corretta solo con pesi >= 0).
    """
    dist = [INF] * len(adj)
    dist[s] = 0
    heap = [(0, s)]
    chiusi = set()
    push, pop, rilassati = 1, 0, 0
    while heap:
        d, u = heapq.heappop(heap)
        pop += 1
        if definitivo:
            if u in chiusi:
                continue
            chiusi.add(u)
        elif d > dist[u]:  # voce vecchia: u e' gia' migliorato
            continue
        for v, w in adj[u]:
            if d + w < dist[v]:
                dist[v] = d + w
                rilassati += 1
                heapq.heappush(heap, (dist[v], v))
                push += 1
    return dist, {"push": push, "pop": pop, "rilassamenti": rilassati}


def bellman_ford(n, archi, s=None):
    """archi = lista di (u, v, peso). Se s e' None parte da 0 su ogni vertice
    (trova un ciclo negativo ovunque). Ritorna dist, ciclo_negativo, passate, controlli."""
    dist = [INF] * n if s is not None else [0] * n
    if s is not None:
        dist[s] = 0
    controlli = 0
    for passata in range(1, n + 1):
        cambiato = False
        for u, v, w in archi:
            controlli += 1
            if dist[u] + w < dist[v]:
                dist[v] = dist[u] + w
                cambiato = True
        if not cambiato:
            return dist, False, passata, controlli
    return dist, True, n, controlli  # ancora cambi alla passata n: ciclo negativo


def floyd_warshall(n, peso):
    """peso[i][j] = peso dell'arco o INF; diagonale 0. Ritorna matrice e passi."""
    d = [riga[:] for riga in peso]
    passi = 0
    for k in range(n):
        dk = d[k]
        for i in range(n):
            dik = d[i][k]
            di = d[i]
            for j in range(n):
                passi += 1
                if dik + dk[j] < di[j]:
                    di[j] = dik + dk[j]
    return d, passi

import random
import networkx as nx
from sp import dijkstra, INF

rng = random.Random(31)
for prova in range(3000):
    n = rng.randint(1, 30)
    G = nx.DiGraph()
    G.add_nodes_from(range(n))
    for _ in range(rng.randint(0, n * 3)):
        u, v = rng.randrange(n), rng.randrange(n)
        if u != v:
            G.add_edge(u, v, weight=rng.randint(0, 20))
    adj = [[(v, G[u][v]["weight"]) for v in G[u]] for u in range(n)]
    s = rng.randrange(n)
    ref = nx.single_source_dijkstra_path_length(G, s)
    for modo in (True, False):
        dist, _ = dijkstra(adj, s, definitivo=modo)
        assert {v: d for v, d in enumerate(dist) if d < INF} == ref
print("3000 digrafi con pesi 0..20: entrambe le varianti uguali a networkx", nx.__version__)

print("V       E        push     pop      rilassamenti  push <= E+1")
for n in (1000, 10000, 100000):
    e = 4 * n
    adj = [[] for _ in range(n)]
    for _ in range(e):
        u, v = rng.randrange(n), rng.randrange(n)
        adj[u].append((v, rng.randint(1, 100)))
    _, c = dijkstra(adj, 0)
    print(f"{n:<7} {e:<8} {c['push']:<8} {c['pop']:<8} {c['rilassamenti']:<13} {c['push'] <= e + 1}")

print()
print("caso peggiore: k vertici u_i, ognuno collegato a tutti gli m vertici v")
print("k=m   E        push     push/E")
for k in (30, 60, 120):
    n = 1 + k + k  # s, u_1..u_k, v_1..v_k
    adj = [[] for _ in range(n)]
    for i in range(1, k + 1):
        adj[0].append((i, i))  # dist(u_i) = i
        for j in range(k + 1, n):
            adj[i].append((j, 2 * k - 2 * i))  # dist(u_i) + peso = 2k - i: cala con i
    e = sum(len(x) for x in adj)
    dist, c = dijkstra(adj, 0)
    print(f"{k:<5} {e:<8} {c['push']:<8} {c['push'] / e:.3f}")

# Output:
# 3000 digrafi con pesi 0..20: entrambe le varianti uguali a networkx 3.7
# V       E        push     pop      rilassamenti  push <= E+1
# 1000    4000     1268     1268     1267          True
# 10000   40000    12878    12878    12877         True
# 100000  400000   129181   129181   129180        True
#
# caso peggiore: k vertici u_i, ognuno collegato a tutti gli m vertici v
# k=m   E        push     push/E
# 30    930      931      1.001
# 60    3660     3661     1.000
# 120   14520    14521    1.000

Entrambe le varianti coincidono con il riferimento su tutti e 3.000 i casi. Nei grafi casuali da 1.000, 10.000 e 100.000 vertici (con quattro archi per vertice) gli inserimenti nello heap sono circa 1,27 per vertice (1.268, 12.878, 129.181): si inserisce solo quando una stima migliora, e in un grafo casuale poche stime migliorano più di una volta. Il limite teorico è E + 1, perché ogni arco può migliorare una stima al massimo una volta (all'estrazione del suo estremo di partenza), e la riga push <= E+1 lo conferma su queste istanze.

Il caso peggiore si costruisce: k vertici u_i con distanza i, ciascuno collegato a ognuno di k vertici v con peso 2k - 2i, così che ogni nuovo u_i estratto migliora tutti i v. Qui gli inserimenti raggiungono E + 1 (931 su 930 archi per k = 30, 14.521 su 14.520 per k = 120): ogni arco costa un inserimento. Ogni operazione sullo heap costa O(log di quanti elementi ci sono), e gli elementi sono al più E + 1, quindi il costo totale è O(E log V): la pagina sui cammini minimi di Algorithms, 4th edition dichiara per Dijkstra con pesi non negativi un tempo proporzionale a E log V nel caso peggiore, e il conteggio sopra ne mostra l'origine, un inserimento per arco.

Dove Dijkstra sbaglia: un arco negativo

L'argomento di Dijkstra («un percorso più lungo non può che costare di più») cade appena un arco può costare meno di zero. Il controesempio minimo: da s si va ad a con costo 1 e a b con costo 4; da b ad a l'arco costa -5, quindi il vero costo di a è 4 - 5 = -1; da a si va a c con costo 1, quindi il vero costo di c è 0. Dijkstra chiude a con costo 1 prima di esaminare b, e quando scopre il percorso migliore non riapre a: il costo di c resta 2.

import random
from sp import dijkstra, bellman_ford, INF

# s -> a costa 1, s -> b costa 4, b -> a costa -5, a -> c costa 1:
# il vero costo di a e' 4 - 5 = -1, quello di c e' 0
archi = [(0, 1, 1), (0, 2, 4), (2, 1, -5), (1, 3, 1)]
adj = [[] for _ in range(4)]
for u, v, w in archi:
    adj[u].append((v, w))
print("Dijkstra classico  :", dijkstra(adj, 0, definitivo=True)[0])
print("Dijkstra con ricontrollo:", dijkstra(adj, 0, definitivo=False)[0])
print("Bellman-Ford       :", bellman_ford(4, archi, 0)[0])

# quanto spesso sbaglia? grafi aciclici (quindi senza cicli negativi), pesi da -5 a 10
rng = random.Random(8)
errori = {True: 0, False: 0}
for prova in range(3000):
    n = rng.randint(2, 12)
    archi = [(u, v, rng.randint(-5, 10)) for u in range(n) for v in range(u + 1, n) if rng.random() < 0.4]
    adj = [[] for _ in range(n)]
    for u, v, w in archi:
        adj[u].append((v, w))
    ref = bellman_ford(n, archi, 0)[0]
    for modo in (True, False):
        if dijkstra(adj, 0, definitivo=modo)[0] != ref:
            errori[modo] += 1
print("su 3000 grafi aciclici con pesi negativi, risposte sbagliate:")
print("  classico:", errori[True], " con ricontrollo:", errori[False])

# Output:
# Dijkstra classico  : [0, -1, 4, 2]
# Dijkstra con ricontrollo: [0, -1, 4, 0]
# Bellman-Ford       : [0, -1, 4, 0]
# su 3000 grafi aciclici con pesi negativi, risposte sbagliate:
#   classico: 105  con ricontrollo: 0

La versione classica dà c = 2 contro il valore vero 0; la variante con ricontrollo e Bellman-Ford danno 0. Il secondo numero dice quanto spesso succede: su 3.000 grafi aciclici casuali (aciclici, così nessun ciclo negativo rende il problema privo di soluzione) con pesi da -5 a 10, la versione classica sbaglia 105 volte, quella con ricontrollo mai. Questo non rende la variante con ricontrollo una soluzione: la pagina di Sedgewick e Wayne vincola la garanzia di E log V alla regola per cui ogni vertice entra una sola volta, e senza quella regola non ho misurato un limite al suo caso peggiore. La scelta onesta con pesi negativi è un altro algoritmo.

Bellman-Ford: rilassare tutti gli archi, V - 1 volte

Bellman-Ford non sceglie un ordine intelligente: rilassa tutti gli archi, in qualunque ordine, e ripete. Dopo la prima passata le stime di tutti i vertici raggiungibili con un solo arco sono giuste; dopo la seconda, quelle raggiungibili con due; un percorso minimo senza cicli ha al massimo V - 1 archi, quindi dopo V - 1 passate tutto è giusto. Se alla V-esima passata una stima cambia ancora, esiste un ciclo di costo negativo raggiungibile, e le distanze non hanno senso perché girare nel ciclo le abbassa all'infinito. La funzione di sp.py si ferma presto se una passata non cambia nulla, e restituisce passate e controlli fatti.

Il listato prova due cose. Su 3.000 digrafi casuali con pesi da -4 a 9 confronta il riconoscimento del ciclo negativo con networkx.negative_edge_cycle (partendo da 0 su ogni vertice, per trovare cicli ovunque) e, dove non c'è ciclo, le distanze con single_source_bellman_ford_path_length. Poi misura il costo su una catena di V vertici, con gli archi elencati in ordine peggiore (dall'ultimo al primo, così ogni passata sistema un vertice solo) e in ordine casuale.

import random
import networkx as nx
from sp import bellman_ford, INF

rng = random.Random(2)
con_ciclo = 0
for prova in range(3000):
    n = rng.randint(2, 14)
    G = nx.DiGraph()
    G.add_nodes_from(range(n))
    for _ in range(rng.randint(0, n * 3)):
        u, v = rng.randrange(n), rng.randrange(n)
        if u != v:
            G.add_edge(u, v, weight=rng.randint(-4, 9))
    archi = [(u, v, G[u][v]["weight"]) for u, v in G.edges]
    _, ciclo, _, _ = bellman_ford(n, archi)  # partenza da ogni vertice
    assert ciclo == nx.negative_edge_cycle(G)
    if ciclo:
        con_ciclo += 1
        continue
    s = rng.randrange(n)
    dist, _, _, _ = bellman_ford(n, archi, s)
    ref = nx.single_source_bellman_ford_path_length(G, s)
    assert {v: d for v, d in enumerate(dist) if d < INF} == ref
print("3000 digrafi con pesi -4..9:", con_ciclo, "con ciclo negativo, riconosciuti come da networkx", nx.__version__, "; negli altri le distanze coincidono")

print("catena di V vertici, archi elencati nell'ordine peggiore o in ordine casuale")
print("V     E      passate(peggiore)  controlli    V*E       passate(casuale)")
for n in (50, 100, 200):
    catena = [(i, i + 1, 1) for i in range(n - 1)]
    _, _, p1, c1 = bellman_ford(n, catena[::-1], 0)
    mescolata = catena[:]
    rng.shuffle(mescolata)
    _, _, p2, _ = bellman_ford(n, mescolata, 0)
    print(f"{n:<5} {n - 1:<6} {p1:<18} {c1:<12} {n * (n - 1):<9} {p2}")

# Output:
# 3000 digrafi con pesi -4..9: 682 con ciclo negativo, riconosciuti come da networkx 3.7 ; negli altri le distanze coincidono
# catena di V vertici, archi elencati nell'ordine peggiore o in ordine casuale
# V     E      passate(peggiore)  controlli    V*E       passate(casuale)
# 50    49     50                 2450         2450      28
# 100   99     100                9900         9900      51
# 200   199    200                39800        39800     101

682 dei 3.000 grafi hanno un ciclo negativo e vengono tutti riconosciuti come fa networkx; negli altri le distanze coincidono. Nella catena con ordine peggiore le passate sono V (le V - 1 necessarie più quella che non cambia nulla) e i controlli sono esattamente V per E: 2.450 per V = 50, 9.900 per V = 100, 39.800 per V = 200. È il caso che dà O(VE). Con gli stessi archi in ordine casuale le passate scendono a circa metà di V (28, 51, 101), perché una passata propaga lontano quando l'ordine è favorevole; resta una dipendenza lineare da V, e nessuna garanzia migliore del caso peggiore.

Floyd-Warshall: tutte le coppie con tre cicli

Se servono le distanze fra tutte le coppie di vertici, ripetere un algoritmo da ogni sorgente funziona, ma Floyd-Warshall (Floyd, Algorithm 97: Shortest path, 1962, sulla base del teorema di Warshall sulle matrici booleane, 1962) usa un'idea più compatta: considera i vertici uno alla volta come possibili tappe intermedie. Dopo aver trattato i vertici da 0 a k, d[i][j] è il costo minimo da i a j usando come tappe solo vertici fra 0 e k. Per passare da k - 1 a k basta una domanda per ogni coppia: conviene passare da k? Cioè d[i][k] + d[k][j] < d[i][j]. L'ordine dei cicli conta: k deve essere il più esterno, altrimenti la tappa k non è ancora «disponibile» quando serve.

La matrice di partenza ha 0 in diagonale, il peso dove c'è un arco e infinito altrove. Se alla fine un elemento della diagonale è negativo, esiste un ciclo negativo: il vertice torna a sé stesso a costo minore di zero. Il listato verifica entrambe le cose contro networkx (matrici uguali a floyd_warshall, ciclo riconosciuto come negative_edge_cycle) su 1.500 digrafi con pesi da -3 a 9, e conta i passi del ciclo più interno su grafi da 20, 40 e 80 vertici.

import random
import networkx as nx
from sp import floyd_warshall, INF

def matrice(n, G):
    m = [[INF] * n for _ in range(n)]
    for i in range(n):
        m[i][i] = 0
    for u, v, w in G.edges(data="weight"):
        m[u][v] = min(m[u][v], w)
    return m

rng = random.Random(4)
ciclo_neg = 0
for prova in range(1500):
    n = rng.randint(1, 12)
    G = nx.DiGraph()
    G.add_nodes_from(range(n))
    for _ in range(rng.randint(0, n * 3)):
        u, v = rng.randrange(n), rng.randrange(n)
        if u != v:
            G.add_edge(u, v, weight=rng.randint(-3, 9))
    d, _ = floyd_warshall(n, matrice(n, G))
    negativo = any(d[i][i] < 0 for i in range(n))  # un vertice che torna da se' a costo < 0
    assert negativo == nx.negative_edge_cycle(G)
    if negativo:
        ciclo_neg += 1
        continue
    ref = dict(nx.floyd_warshall(G))
    assert all(d[i][j] == ref[i][j] for i in range(n) for j in range(n))
print("1500 digrafi: matrici uguali a networkx", nx.__version__, "e", ciclo_neg, "cicli negativi riconosciuti dalla diagonale")

print("V     passi del ciclo interno   V^3")
for n in (20, 40, 80):
    G = nx.gnp_random_graph(n, 0.2, seed=3, directed=True)
    nx.set_edge_attributes(G, 1, "weight")
    _, passi = floyd_warshall(n, matrice(n, G))
    print(f"{n:<5} {passi:<25} {n ** 3}")

# Output:
# 1500 digrafi: matrici uguali a networkx 3.7 e 207 cicli negativi riconosciuti dalla diagonale
# V     passi del ciclo interno   V^3
# 20    8000                      8000
# 40    64000                     64000
# 80    512000                    512000

I passi sono esattamente V al cubo: 8.000, 64.000 e 512.000. Il costo non dipende da quanti archi ci sono, perché il ciclo triplo scorre ogni coppia per ogni tappa: su un grafo con pochi archi Floyd-Warshall costa lo stesso che su uno denso, ed è il motivo per cui conviene solo quando V è piccolo o il grafo è denso. In memoria occupa V al quadrato celle, la stessa matrice di adiacenza dell'articolo 6.

Quale scegliere

Le tre complessità misurate stanno in una tabella; i costi sono quelli dei conteggi sopra, non tempi.

La lettura corretta della tabella parte dalla domanda sui pesi: se sono tutti non negativi (distanze, tempi, costi reali) Dijkstra è la scelta, perché è l'unico dei tre il cui costo dipende da E log V invece che da V per E. Se possono essere negativi (un guadagno in una catena di scambi, uno sconto), serve Bellman-Ford per una sorgente, e il suo controllo del ciclo negativo ti dice se la domanda ha senso. Se servono tutte le coppie e V è piccolo, Floyd-Warshall è il più breve da scrivere e da controllare.

Non ho implementato né misurato A*, le gerarchie di contrazione o la ricerca bidirezionale: li nomino perché esistono varianti più veloci di Dijkstra per grafi stradali molto grandi, ma ogni affermazione su quanto siano più veloci richiederebbe una fonte e una misura che qui non ci sono.

Che cosa sai fare adesso

Sai rilassare un arco, sai scrivere Dijkstra con heapq senza l'operazione «abbassa la priorità», e sai costruire il grafo di quattro vertici in cui sbaglia perché un peso è negativo. Sai usare Bellman-Ford per riconoscere un ciclo negativo, e Floyd-Warshall per avere tutte le distanze con un ciclo triplo. Sai dichiarare ognuna delle tre complessità con un conteggio che il lettore può rifare. Il prossimo e ultimo articolo, sulla programmazione dinamica, riprende l'idea di Floyd-Warshall (una tabella riempita in un ordine che garantisce che le risposte parziali ci siano già) e la applica a problemi che non sono grafi. Per tornare indietro, l'articolo 5 è quello sugli alberi.

Nella stessa serie trovi anche la complessità e Big-O, gli algoritmi di ordinamento, la ricerca lineare, binaria e con tabelle hash e liste, pile, code e deque.