LabASD - Algoritmo di Bellman-Ford

Implementare l’algoritmo di Bellman-Ford per il calcolo del cammini minimi da singola sorgente in un grafo orientato pesato. L’algoritmo richiede un nodo di partenza \(s\), e al termine calcola la distanza \(d[v]\) di ogni nodo \(v\) dalla sorgente, \(v \in \{0, \ldots, n-1\}\). Se un nodo non è raggiungibile, si avrà \(d[v] = +\infty\). Se il grafo contiene un ciclo negativo che è raggiungibile dalla sorgente, il programma si deve interrompere stampando un opportuno messaggio di errore.

L’algoritmo di Bellman-Ford opera eseguendo dei passi di rilassamento. Dato un arco orientato \((u, v)\) di peso \(w\), un passo di rilassamento consiste nell’operazione seguente:

if (d[u] + w < d[v]) {
    d[v] = d[u] + w;
    p[v] = u;
}

dove p[v] indica il predecessore del nodo \(v\) lungo il cammino minimo dalla sorgente a \(v\); se \(v\) non ha predecessore (ad es., perché è il nodo sorgente o perché non è raggiungibile dalla sorgente) si deve porre p[v] = -1 oppure p[v] = NODE_UNDEF. L’array p[] serve per poter ricostruire i cammini minimi.

L’operazione di rilassamento riduce la stima \(d[v]\) della distanza del nodo \(v\) dalla sorgente, nel caso in cui si scopra che attraversando l’arco \((u, v)\) si raggiunga \(v\) con una distanza inferiore. Le stime iniziali delle distanze sono tutte \(+\infty\), ad esclusione della sorgente \(s\) per cui si ha d[s] = 0.

Per rappresentare \(+\infty\) si può usare il simbolo HUGE_VAL che è associato ad un valore di tipo float o double definito in <math.h>¹. Questo simbolo è parte dello standard C90, per cui deve essere presente in tutti i compilatori conformi.

Questo programma memorizza l’albero dei cammini minimi in due modi: come array di predecessori p[v] di lunghezza \(n\), e in un array di puntatori ad archi sp[] sempre di lunghezza \(n\). Specificamente, p[v] è un intero e rappresenta l’identificativo (indice) del nodo che precede \(v\) lungo il cammino minimo dalla sorgente a \(v\); se \(v\) non è raggiungibile dalla sorgente, oppure coincide con la sorgente e quindi non ha predecessore, si pone p[v] = -1 oppure p[v] = NODE_UNDEF (il simbolo NODE_UNDEF vale -1 ed è definito all’inizio del programma). sp[v] è un puntatore all’arco entrante nel nodo \(v\) lungo il cammino minimo dalla sorgente a \(v\); se \(v\) non è raggiungibile dalla sorgente, oppure coincide con la sorgente, si pone sp[v] = NULL; si faccia riferimento alla Figura 1.

Poiché l’array sp[] non viene mai sfruttato da questo programma, suggerisco di iniziare lasciandolo indefinito.

Il programma riceve da uno a tre parametri sulla riga di comando:

• il primo parametro è il nome del file di input; il formato dell’input è quello utilizzato dal tipo grafo;

• il secondo parametro è un intero che indica il nodo sorgente (se manca, si assume 0 come sorgente);

• il terzo parametro è un intero che indica il nodo destinazione. Si noti che il programma calcola sempre le distanze di tutti i nodi dalla sorgente, per cui il nodo destinazione influenza solo ciò che viene stampato a video ma non il costo asintotico dell’algoritmo. Se la destinazione non viene specificata, stampa le distanze e i cammini minimi verso tutti i nodi del grafo.

Per compilare:

    gcc -std=c90 -Wall -Wpedantic graph.c bellman-ford.c -o bellman-ford

Per stampare le distanze e i cammini minimi dal nodo 0:

    ./bellman-ford graph10.in 0

su Windows:

    .\bellman-ford graph10.in 0

Per stampare la distanza e il cammino minimo dal nodo 3 al nodo 4:

    ./bellman-ford graph10.in 3 4

su Windows:

    .\bellman-ford graph10.in 3 4

Attenzione ad un bug di alcuni processori/compilatori

La procedura di rilassamento può causare problemi su alcune combinazioni di processori/compilatori. Consideriamo un arco \((u, v)\) di peso \(w\). Ci aspetteremmo che nel frammento di codice seguente l’asserzione risulti sempre verificata:

if (d[u] + w < d[v]) {
   d[v] = d[u] + w;
   assert(d[v] == d[u] + w);
}

Purtroppo, in certi casi rari può capitare che l’asserzione venga violata! Questo comporta che l’algoritmo non converga mai anche se non sono presenti cicli negativi, perché continuamente rilassati gli stessi archi anche quando si è determinata la distanza minima della destinazione. Il risultato è che il programma riporta la presenza di cicli negativi anche quando non ce ne sono.

Ho riscontrato il problema su un processore Intel Atom con gcc 7.x/8.x oppure con tcc (un altro compilatore per sistemi Unix), ma non con clang (il compilatore usato di default su MacOSX). Curiosamente, sul processore Atom non funziona nemmeno una implementazione Pascal dell’algoritmo, compilata con Free Pascal 3.x. Ciò fa sospettare che il problema sia legato in qualche modo all’hardware.

L’errore si è manifestato anche durante il laboratorio su un portatile dotato di processore Intel i7 usando una vecchia versione di gcc (4.x). Il problema non sembra manifestarsi con processori recenti usando versioni recenti di gcc (ho provato: Xeon, i7, e ARMv7l su una scheda Raspberry Pi). Riassumendo:

Il problema sembra quello descritto in https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323 e potrebbe derivare dal fatto che il coprocessore matematico 8087 non era totalmente conforme allo standard IEEE 754; in particolare, i registri in virgola mobile avevano (e hanno tutt’ora) una ampiezza maggiore di 32/64 bit, per cui la somma viene calcolata con un numero maggiore di cifre, mentre viene memorizzata con meno (excess precision).

Possibili soluzioni:

• Se si usa gcc è possibile forzare l’uso delle istruzioni SSE compilando con i flag -msse2 -mfpmath=sse. In base a quanto scritto in https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323#c109 le istruzioni SSE sono conformi IEEE 754, e il risultato della computazione dovrebbe essere corretto.

• Riscrivere la procedura di rilassamento introducendo una variabile temporanea (questa è probabilmente la soluzione più semplice):

  const double dnew = d[u] + w;
  if (dnew < d[v]) {
    d[v] = dnew;
  }

Sul processore Atom, clang (che produce codice funzionante) sembra generare istruzioni SIMD per il calcolo in virgola mobile, che non vengono usate di default da gcc.

Nella maggior parte dei processori recenti con una versione recente di gcc non dovrebbe servire nessuna delle “pezze” di cui sopra.

Suggerimento per una possibile ottimizzazione

Nella versione base l’algoritmo di Bellman-Ford effettua \((n-1)\) fasi di rilassamento; in ogni fase si esaminano tutti gli archi del grafo, per cui il costo dell’algoritmo è \(\Theta(nm)\).

Conviene però adottare una semplice ottimizzazione: se al termine di una fase di rilassamento le distanze dei nodi non cambiano, allora l’algoritmo può terminare immediatamente perché i rilassamenti successivi non modificheranno le distanze.

Questo comportamento si può realizzare con una modifica alla versione base dell’algoritmo vista a lezione. Si noti che la funzione relax() proposta in questo programma restituisce 1 se la distanza del nodo destinazione viene modificata, 0 altrimenti. Ciò è stato fatto proprio per consentire l’ottimizzazione di cui stiamo parlando.

int i, v, updated;

/* Inizializzazione non riportata... */

for (updated=1, i=0; i<n && updated; i++) {
  updated = 0;
  for (v=0; v<n; v++) {
    const Edge *e;
    for (e = graph_adj(g, v); e != NULL; e = e->next) {
      updated |= relax(e, d, p, sp);
    }
  }
}

Il codice precedente effettua un passo di rilassamento in più: nello pseudocodice il ciclo esterno viene eseguito \((n-1)\) volte, non \(n\) come nel frammento sopra. Questo è voluto, e serve per verificare se ci sono cicli negativi. Nei lucidi mostrati a lezione viene fatto un ciclo in più per controllare la presenza di cicli negativi. Nel codice sopra, invece, basta controllare se alla fine si ha updated == 1. Se così fosse, significherebbe che sono state fatte \(n\) iterazioni, ed al termine dell’ultima sono state modificate delle distanze.

Con l’ottimizzazione di cui sopra è possibile gestire grafi più grandi, come quelli nella Tabella 1 che rappresentano porzioni delle reti stradali di alcuni stati americani e della città di Roma. Viste le dimensioni degli input, conviene disattivare la stampa dei cammini minimi (commentare la chiamata alla funzione print_path() dentro print_dist()).

I dati sono forniti nel 9th DIMACS implementation challenge–Shortest Paths. Prestare attenzione che i grafi non sono sempre (fortemente) connessi, quindi non tutti i nodi sono raggiungibili dagli tutti altri.

File

• bellman-ford.c
• graph.c
• graph.h
• graph-negative-cycle.in: questo grafo ha cicli negativi, per cui usando il nodo 0 come sorgente il programma deve stampare un messaggio di errore.
• graph10.in (risultato atteso usando il nodo 0 come sorgente)
• graph100.in
• graph1000.in