Recepty z programátorské kuchařky
Haldy, heapsort a Dijkstrův algoritmus

Dnešní povídání o algoritmech a datových strukturách se bude zabývat jedním z nejznámějších algoritmů: Dijkstrovým algoritmem pro hledání nejkratších cest v grafech. K tomu (a nejenom k tomu) se nám bude hodit šikovná datová strukturka zvaná halda, tak si předvedeme nejdříve ji.

Halda

Halda je datová struktura pro uchovávání množiny čísel (či jakýchkoliv jiných prvků, na kterých máme definováno uspořádání, tj. umíme pro každou dvojici prvků říci, který z nich je menší). Tato datová struktura obvykle podporuje následující operace: přidání nového prvku, nalezení nejmenšího prvku a odebrání nejmenšího prvku. My si ukážeme jednoduchou implementaci haldy, která bude při uložení N prvků potřebovat čas O(log N) na přidání či odebrání jednoho prvku a O(1) (tj. konstantní) na zjištění hodnoty nejmenšího prvku.

Naše implementace bude vypadat následovně: Pokud halda obsahuje N prvků, uložíme její prvky do pole na pozice 1 až N. Prvek na pozici k bude mít dva následníky, a to prvky na pozicích 2k a 2k+1; samozřejmě, pokud je k velké, a tedy např. 2k+1>N, má takový prvek jen jednoho či dokonce žádného následníka. Naopak prvek na pozici ⌊k/2⌋ nazveme předchůdcem prvku na pozici k. Ti z vás, kteří znají binární stromy, v tomto jistě rozpoznali způsob, jak v poli uchovávat úplné binární stromy (následníci jsou synové a předchůdci otcové v obvyklé stromové terminologii, prvek č. 1 je kořen stromu).

Prvky haldy však v poli neuchováváme v úplně libovolném pořadí. Chceme, aby platilo, že každý prvek je menší nebo roven všem svým následníkům. Naše halda tedy může vypadat např. takto:

Tomu odpovídá tento strom:

Z toho, co jsme si právě popsali, je jasné, že nejmenší prvek je uložen na pozici s indexem 1, a tedy můžeme snadno v konstantním čase zjistit jeho hodnotu. Ještě prozradíme, jak lze prvky do haldy rychle přidávat a odebírat:

Jestliže halda obsahuje N prvků, pak nový prvek, říkejme mu třeba x, přidáme na konec pole, tj. na pozici s indexem N+1. Nyní x porovnáme s jeho předchůdcem. Pokud je jeho předchůdce menší, je vše v pořádku a jsme hotovi. V opačném případě x s jeho předchůdcem prohodíme. Tím jsme problém napravili, ale nyní může být x menší než jeho nový předchůdce. To lze napravit dalším prohozením a tak budeme pokračovat dále, než se buďto dostaneme do situace, kdy už je x větší nebo rovno svému předchůdci, nebo „vybublá“ až do kořene haldy, kde už žádného předchůdce nemá. Protože se v každém kroku pozice, na níž se prvek x právě nachází, zmenší alespoň na polovinu, provedeme dohromady nejvýše O(log N) výměn, a tedy spotřebujeme čas O(log N).

Odebírání nejmenšího prvku probíhá podobně: Prvek z poslední pozice (tj. z pozice N) přesuneme na pozici 1, tedy místo minima. Místo s předchůdci jej však porovnáme s jeho následníky, a v případě, že je větší než některý z jeho následníků, opět je prohodíme (pokud je větší než oba následníci, prohodíme ho s menším z nich). A protože se nám v každém kroku index „bublajícího“ prvku v poli alespoň zdvojnásobí, opět spotřebujeme čas O(log N).

Jako cvičení si rozmyslete, že v čase O(log N) lze z haldy smazat dokonce libovolný prvek, pokud si ovšem pamatujeme, kde se v haldě nachází. Také můžeme prvek ponechat a jen změnit jeho hodnotu.

Ještě si předvedeme program:

halda = []

def nejmensi():
  return halda[0]

def swap(a,b):
  (halda[a], halda[b]) = (halda[b], halda[a])

def vloz(prvek):
  # Vložíme prvek na konec
  halda.append(prvek)
  i = len(halda) - 1
  while (i > 0) and (halda[i//2] > halda[i]):
    swap(i, i//2)
    i = i//2

def smaz_nejmensi():
  N = len(halda)
  nejmensi = halda[0]
  halda[0] = halda[N - 1]
  # Odstraníme prvek z konce
  halda.pop()
  N -= 1
  i = 0
  while 2*i < N:
    j = i
    if halda[j] > halda[2*i]:
      j = 2*i
    if (2*i+1 < N) and (halda[j]>halda[2*i+1]):
      j = 2*i+1
    if i == j:
      break
    swap(i, j)
    i = j
  return nejmensi

var halda: array[1..MAX] of integer;
    N: integer; { počet prvků v haldě }

function nejmensi: integer;
begin
  nejmensi:=halda[1]
end;

procedure vloz(prvek: integer);
var i, x: integer;
begin
  N:=N+1; i:=N;
  halda[i]:=prvek;
  while (i>1) and (halda[i div 2]>halda[i])
  do begin
    x:=halda[i div 2];
    halda[i div 2]:=halda[i];
    halda[i]:=x;
    i:=i div 2
  end
end;

procedure smaz_nejmensi;
var i, j, x: integer;
begin
  halda[1]:=halda[N];
  N:=N-1; i:=1;
  while 2*i<=N do begin
    j:=i;
    if halda[j]>halda[2*i] then j:=2*i;
    if (2*i+1<=N) and (halda[j]>halda[2*i+1])
       then j:=2*i+1;
    if i=j then break;
    x:=halda[i]; halda[i]:=halda[j];
    halda[j]:=x;
    i:=j
  end
end;

HeapSort

Když už máme k dispozici haldu, můžeme pomocí ní například snadno třídit čísla. Máme-li N čísel, která chceme setřídit, vytvoříme si z nich nejprve haldu o N prvcích (například postupným vkládáním do prázdné haldy), načež z ní budeme postupně N-krát odebírat nejmenší prvek. Tím získáme prvky původního pole v rostoucím pořadí. Celkově provedeme N vložení, N nalezení minima a N smazání. To vše dohromady stihneme v čase O(N log N).

Než si ukážeme program, přidáme ještě dva triky, které nám implementaci značně usnadní. Předně si vše uložíme do jednoho pole – to bude při plnění haldy obsahovat na svém začátku haldu a na konci zbytek vstupního pole, přitom zbytek pole se bude postupně zmenšovat a uvolňovat tak místo haldě; naopak v druhé polovině algoritmu budeme zmenšovat haldu a do volného prostoru ukládat setříděné prvky. K tomu se nám bude hodit získávat prvky v opačném pořadí, proto si upravíme haldu tak, aby udržovala nikoliv minimum, nýbrž maximum.

Druhý trik spočívá v tom, že nebudeme haldu vytvářet postupným vkládáním, nýbrž naopak zabubláváním prvků (podobným, jako děláme při mazání minima) od konce. Všimněte si, že takto také získáme správné nerovnosti mezi prvky a jejich následníky, a dokonce tak zvládneme celou haldu vytvořit v lineárním čase (proč to tak je, si zkuste dokázat sami, stačí si uvědomit, kolikrát zabubláváme které prvky). Zbytek třídění bohužel nadále zůstává O(N log N).

Tomuto algoritmu se obvykle říká HeapSort (čili třídění haldou) a je jedním z mála známých rychlých třídících algoritmů, které nepotřebují pomocnou paměť.

# Zabublání prvku dolů:
#   N určuje část pole vyhrazenou haldě
#   i je index zabublávaného prvku
def bublej(pole, N, i):
  while 2*i < N:
    j = 2*i
    if (j+1 < N) and (pole[j+1] > pole[j]):
      j = j+1
    if pole[i] >= pole[j]:
      break
    (pole[i], pole[j]) = (pole[j], pole[i])
    i = j

def heapsort(pole):
  N = len(pole)
  for i in range(N//2, -1, -1):
    bublej(pole, N, i)
  # Výběr maxima a jeho přesun nakonec
  for i in range(N - 1, 0, -1):
    (pole[0], pole[i]) = (pole[i], pole[0])
    # Dál už bubláme jen na zbytku pole
    bublej(pole, i, 0)
  return pole

type Pole = array[1..MAXN] of Integer;

procedure HeapSort(var A: Pole);
var i, x: integer;
   procedure bublej(m, i: integer);
   { zabublání prvku: m je velikost haldy,
   i je index zabublávaného prvku }
   var j, x: integer;
   begin
      while 2*i<=m do begin
         j:=2*i;
         if (j<m) and (A[j+1]>A[j]) then j:=j+1;
         if A[i]>=A[j] then break;
         x:=A[i]; A[i]:=A[j]; A[j]:=x;
         i:=j;
      end;
   end;
begin
   for i:=N div 2 downto 1 do bublej(N,i);
   { vybírej maximum }
   for i:=N downto 2 do begin
      x:=A[1]; A[1]:=A[i]; A[i]:=x;
      bublej(i-1, 1);
   end;
end;

Dijkstrův algoritmus

Nyní již konečně k slíbenému Dijkstrovu algoritmu. Tento algoritmus dostane orientovaný graf s hranami ohodnocenými nezápornými čísly (viz kuchařka o grafech) a nalezne v něm nejkratší cestu mezi dvěma zadanými vrcholy. Ve skutečnosti tento algoritmus dělá o malinko více: Najde totiž nejkratší cesty z jednoho zadaného vrcholu do všech ostatních.

Nechť v₀ je vrchol grafu, ze kterého chceme určit délky nejkratších cest. Budeme si udržovat pole délek zatím nalezených cest z vrcholu v₀ do všech ostatních vrcholů grafu. Navíc u některých vrcholů budeme mít poznamenáno, že cesta nalezená do nich je už ta nejkratší možná. Takovým vrcholům budeme říkat definitivní. Na začátku inicializujeme v poli všechny hodnoty na ∞ kromě hodnoty odpovídající vrcholu v₀, kterou inicializujeme na 0 (délka nejkratší cesty z v₀ do v₀ je 0). V každém kroku algoritmu pak provedeme následující: Vybereme vrchol w, který ještě není definitivní, a mezi všemi takovými vrcholy je délka zatím nalezené cesty do něj nejkratší možná. Vrchol w prohlásíme za definitivní. Dále otestujeme, zda pro nějaký vrchol v cesta z vrcholu v₀ do w a pak po hraně z w do v není kratší, než zatím nalezená cesta z v₀ do v, a je-li tomu tak, upravíme délku zatím nalezené cesty do v. Toto provedeme pro všechny takové vrcholy v. Celý algoritmus skončí, pokud jsou už všechny vrcholy definitivní nebo všechny vrcholy, které nejsou definitivní, mají délku cesty rovnou ∞ (v takovém případě se graf skládá z více nesouvislých částí).

Předtím, než dokážeme, že právě představený algoritmus opravdu nalezne délky nejkratších cest z vrcholu v₀, se zamysleme nad jeho časovou složitostí.

Pro každý z N vrcholů si délku dosud nalezené cesty uchováme v poli. Celý algoritmus provede nejvýše N kroků, protože v každém kroku nám přibyde jeden definitivní vrchol. Ten vybíráme jako minimum z délky aktuální cesty přes všechny dosud nedefinitivní vrcholy, kterých je O(N). V každém kroku musíme zkontrolovat tolik vrcholů v, kolik hran vede z vrcholu w. Počet takových změn pro všechny kroky dohromady je pak nejvýše O(M), kde M je počet hran vstupního grafu. Z toho vyjde časová složitost O(N²+M), čili O(N²), jelikož M je nejvýše N². Tuto implementaci Dijkstrova algoritmu najdete na konci naší kuchařky.

K uchovávání délek dosud nalezených nejkratších cest můžeme ovšem použít haldu. Ta bude na začátku obsahovat N prvků a v každém kroku se počet jejích prvků sníží o jeden: Nalezneme a smažeme nejmenší prvek, to zvládneme v čase O(log N), a případně upravíme délky nejkratších cest do sousedů právě zpracovávaného vrcholu. To pro každou hranu trvá rovněž O(log N), celkově za všechny hrany tedy O(M log N). Z toho vyjde celková časová složitost algoritmu O((N+M) log N), a to je pro „řídké“ grafy (tedy grafy s M≪N²) výrazně lepší.

Vraťme se nyní k důkazu správnosti Dijkstrova algoritmu. Ukážeme, že po každém kroku algoritmu platí následující tvrzení: Nechť A je množina definitivních vrcholů. Pak délka dosud nalezené cesty z v₀ do v (v je libovolný vrchol grafu) je délka nejkratší cesty v₀v₁… v_k v takové, že všechny vrcholy v₀,v₁,… ,v_k jsou v množině A. Tvrzení dokážeme indukcí dle počtu kroků algoritmu, které již proběhly. Tvrzení zřejmě platí před a po prvním kroku algoritmu. Nechť w je vrchol, který byl v předchozím kroku prohlášen za definitivní. Uvažme nejprve nějaký vrchol v, který je definitivní. Pokud v=w, tvrzení je triviální. V opačném případě ukážeme, že existuje nejkratší cesta z v₀ do v přes vrcholy z A, která nepoužívá vrchol w. Označme D délku cesty z v₀ do v přes vrcholy A bez vrcholu w. Protože v každém kroku vybíráme vrchol s nejmenším ohodnocením a ohodnocení vybraných vrcholů v jednotlivých krocích tvoří neklesající posloupnost (váhy hran jsou nezáporné!), tak délka cesty z v₀ do w přes vrcholy z A je alespoň D. Ale potom délka libovolné cesty z v₀ do v přes w používající vrcholy z A je alespoň D. Z volby D pak víme, že existuje nejkratší cesta z v₀ do v přes vrcholy z A, která nepoužívá vrchol w.

Nyní uvažme takový vrchol v, který není definitivní. Nechť v₀v₁… v_k v je nejkratší cesta z v₀ do v taková, že všechny vrcholy v₀,v₁,… ,v_k jsou v množině A. Pokud v_k=w, pak jsme ohodnocení v změnili na délku této cesty v právě proběhlém kroku. Pokud v_k≠ w, pak v₀v₁,… ,v_k je nejkratší cesta z v₀ do v_k přes vrcholy z množiny A, a tedy můžeme předpokládat, že žádný z vrcholů v₁,… ,v_k není w (podle toho, co jsme si rozmysleli na konci minulého odstavce). Potom se ale délka cesty do v rovnala správné hodnotě už před právě proběhlým krokem.

Vzhledem k tomu, že po posledním kroku množina A obsahuje právě ty vrcholy, do kterých existuje cesta z vrcholu v₀, dokázali jsme, že náš algoritmus funguje správně.

Na závěr ještě poznamenejme, že Dijkstrův algoritmus je možné snadno upravit tak, aby nám kromě určení délky nejkratší cesty i takovou cestu našel: U každého vrcholu si v okamžiku, kdy mu měníme ohodnocení, poznamenáme, ze kterého vrcholu do něj přicházíme. Nejkratší cestu do nějakého vrcholu pak zrekonstruujeme tak, že u posledního vrcholu této cesty zjistíme, který vrchol je předposlední, u předposledního, který je předpředposlední, atd.

Poznámka pro zvídavé: existují i jiné druhy hald, například k-regulární haldy, v nichž má každý prvek k následníků (rozmyslete si, jaká je v takové haldě časová složitost operací a jak nastavit k v závislosti na M a N, aby byl Dijkstrův algoritmus co nejrychlejší), nebo tzv. Fibonacciho halda, která dokáže upravit hodnotu prvku v konstantním čase. S tou pak umíme hledat nejkratší cesty v čase O(M+N log N).

Dan Kráľ, Martin Mareš a Petr Škoda

Implementace Dijkstrova algoritmu s haldou

# Struktura pro vrchol jako dvojici (index, vzdálenost)
# (definuje vlastní porovnávání, aby šlo použít haldu výše)
class vrchol():
  def __init__(self, index, vzdalenost):
    self.index = index
    self.vzdalenost = vzdalenost
  # Předefinování operátoru >
  def __gt__(self, obj):
    return self.vzdalenost.__gt__(obj.vzdalenost)

halda = []
def dijkstra(N, cesty, start):
  final = [False] * N
  vzdalenost = [None] * N
  # Inicializace startu:
  vloz(vrchol(start, 0))
  vzdalenost[start] = 0
  while halda:
    # Vytáhneme z haldy nejmenší nezpracované místo
    v = smaz_nejmensi()
    if final[v.index]:
      continue
    # Označíme místo za použité
    final[v.index] = True
    # Projdeme všechny sousedy
    for (i, delka) in cesty[v.index]:
      if vzdalenost[i] is None or v.vzdalenost + delka < vzdalenost[i]:
        vzdalenost[i] = v.vzdalenost + delka
        vloz(vrchol(i, vzdalenost[i]))
  return vzdalenost