Pokročilé datové struktury a algoritmy
ⓒ Jiří Fišer, MMXX, Přírodovědecká fakulta, UJEP, Ústí nad Labem
1. Asymptotická složitost¶
Časová složitost¶
Při hodnocení algoritmů je klíčovým pojmem asymptotická časová složitost, který umožní klasifikovat algoritmy podle závislosti mezi velikostí úlohy a časem provedení.
Formální zavedení viz například: https://cw.fel.cvut.cz/b182/_media/courses/b6b36dsa/dsa-3-slozitostalgoritmu.pdf.
Přehled běžně se vyskytujích tříd nabízi anglická wikipedia (https://en.wikipedia.org/wiki/Time_complexity).
Pro praxi je třeba znát tato tvrzení a jejich důsledky:
- řád růstu funkcí zohledňuje velké úlohy (velký počet zpracovávaných dat) a nemusí být relevantní pro malé úlohy (sofistikované algoritmy s "lepší" časovou složitostí mohou být na malých datech výrazně pomalejší než algoritmy triviální, například optimálním vyhledávacím algoritmem pro data s počtem položek v řádu jednotek je procházení seznamu a nikoliv hashovací tabulka).
- asymptotická složitost často závisí na vstupních datech. Běžně se uvádí tzv. typická složitost, která je je dosažitelná s vysokou pravděpodobností pro náhodná vstupní data. Pro některé typy vstupních data však může být lepší (například třídění již uspořádaných dat algoritmem insertion sort má složitost $O(n)$ oproti typické $O(n^2)$), ale i horší (některé implementace algoritmu quick sort, který má typicky složitost $O(n\cdot log(n))$, třídí již setříděný seznam se složitostí $O(n^2)$)). Na rozdíl od notace používané v matematické teorie, se v praxi používá jen big O notace.
- běžně sae předpokládá, že všechny elementární operace (porovnání, aritmetické operace, přesuny) mají konstantní časovou složitost (tj. využívají pro ukládání čísel či dalších objektů fixní paměťový prostor. To nemusí platit, pokud se využívá aritmetika s neomezenou přesností nebo řetězce s neomezenou délkou)
Jaká je asympotická složitost následujících operací nad seznamem? (výmaz v tomto případě zahrnuje i posun všech následujících prvků o jednu pozici vlevo)
- výmaz $n-k$ prvku kde $k$ je předem daná konstanta pro níž platí $\forall n, 0 \leq k < n$
- výmaz $n/k$ prvku, kde výmaz $n-k$ prvku kde $k$ je předem daná konstanta pro níž platí $\forall n, 2 \leq k < n$
Paměťová složitost¶
Paměťová složitost nehraje tak klíčovou roli. V praxi, lze rozeznávat následující případy:
- algoritmus vyžaduje konstantní množství paměti bez ohledu na velikost
- algoritmus nepotřebuje kromě vstupu žádnou dodatečnou paměť nebo jen (malou) konstantní či na logaritmu závislou paměť (většina třídících algoritmů, vyhledávání)
- algoritmus potřebuje kopii původních dat resp. pomocnou strukturu s podobnou velikostí
- algoritmus potřebuje ($k$-krát více paměti než vstup, kde $k$ je výrazně větší než 1)
- algorimus využívá pomocnou paměť, jejíž velikost roste rychleji než lineárně
Algoritmy prvních tří skupin jsou aplikovatelné na počítačích, jejichž velikost (operační)paměti je nemusí být řádově větší než velikost vstupních dat. Do této skupiny spadá většina běžně používaných algoritmů a datových struktur.
Algoritmů třetí a skupiny, mohou být navzdory excelentní časové složitosti nepoužitelné pro větší data (a tím být de facto diskvalifikovány), neboť vyžadují odkládání dat na pomalejší vnější paměťová zařízení resp. se nevejdou ani do rozšířené paměti.
Příkladem může být například algoritmus bucket sort, který má časovou složitost $O(n + \frac{n^2}{k} + k)$ a prostorovou $O(n\cdot k)$. Je-li $k$ velké, pak je časová složitost téměř lineární, ale pomocná data se již nemusí do paměti vejít.
Benchmarking¶
Časovou složitost některých algoritmů lze relativně snadno odvodit (například u algoritmů, které procházejí jednoduché kolekce pomocí cyklů), u některých je to již složitější (některé typy rekurzivních algoritmů), a existují však algoritmy, u nichž lze složitost i za použití vyšší matematiky jen přibližně odhadnout.
Pro reálné nasazení je často nutné znát i přesnou závislost času vykonávání na velikosti úlohy. Tuto funkci lze ze zápisu algoritmu odvodit jen velmi přibližně na základě zjištění počtu jednotlivých operací. Reálnou rychlost implementace na konkrétní však silně ovlivňuje také:
- relativní rychlost jednotlivých operací (násobení je na některých platformách stejně rychlé jako sčítání, na jiných [bez harwarové násobičky]) je řádově pomalejší)
- interní paralelismus daný především šířkou zpracovávaných dat, zřetězeným zpracováváním operací (pipelining) a v některých případech i využitím SIMD vektorizace
- velikostí a charakterem jednotlivých úrovní mezipamětí (mezipaměť může výrazně ovlivnit efektiviru některých algoritmů)
Z tohoto důvodu má smysl využívat při implementaci a porovnání algoritmů benchmarking.
Následující kód ukazuje implementuje jednoduchou funkci benchmark, která umožňuje
testovat funkci resp. seznam funkcí pro různé hodnoty $n$ a několik dalších
pomocných funkcí (funkce plot výsledek vizualizuje prostřednictvím,
jednoduchého grafu).
Poznámka: Obsah buňky s kódem je uložen do souboru
algbench.pya může tak být importován jako modul.
%%writefile algbench.py
import collections.abc
from time import perf_counter
def benchmark(algo_function, init, nrange, repetition=100):
if not isinstance(algo_function, collections.abc.Iterable):
algo_function = [algo_function]
nlist = list(nrange)
results = []
for algo in algo_function:
times = []
for n in nlist:
timesum = 0.0
for _ in range(repetition):
input = init(n)
start = perf_counter()
algo(input)
end = perf_counter()
timesum += end-start
times.append(timesum/repetition)
results.append(times)
return (nlist, results)
def qplot(tdata, labels):
import matplotlib.pyplot as plt
for (times, label) in zip(tdata[1], labels):
plt.plot(tdata[0], times, label=label)
plt.xlabel("n")
plt.ylabel("t(s)")
plt.legend()
plt.show()
def nrange(start, n = 10, step = None):
step = step or start
for i in range(start, start + n*step, step):
yield i
Pro otestování provedeme benchmark implementace standardní pythonské funkce max a tří
implementací funkce sum (standardní pythonské přes iterátor a verze z knihovny numpy
přes numpy n-dimenzionální pole, redukce operací plus přes iterátor) a také numpy
implementace funkce sort
from algbench import *
import numpy as np
from functools import reduce
import operator
r = benchmark(
[max, sum, np.sum, lambda input: reduce(operator.add, input, 0.0), np.sort],
lambda n: np.random.uniform(0.0,1.0,n), nrange(10000))
qplot(r, ["max", "sum", "np.sum", "reduce", "sort"])
Výsledek není překvapivý, neboť všechny implementace až na třídění používají
triviální algoritmus, ve kterém
se prochází všechny prvky pole a buď se sčítají, nebo se ukládá průběžné maximum
(rychlejší algoritmus v tomto případě neexistuje).
Časová složitost je tak lineární tj. $O(n)$ (i když směrnice u numpy
verze je tak malá, že se na grafu jeví téměř jako konstantní).
Třídící algoritmus se jeví jako by měl lineární složitost, to je však dáno volbou
rozsahu velikostí polí (skutečbost je samozřejmě jiná viz úkol).
Mírně překvapit může velký rozdíl mezi numpy verzí (implementovanou v jazyce C) a verzemi
vytvořenými přímo v Pythonu (z nichž nejhorší je ta využívající funkcionálního přístupu,
Python není na rozdíl od jazyka Julia nakloněn funkcionálnímu programování).
Pro pole velikostí desetitisíců položek je třídění pomocí numpy funkce stejně rychlé
jako sumace či nalezení maximálního prvku).
Těžko vysvětlitelná je i mírně rychlejší implementace vestavěné funkce max a oproti sum
(v obou se v zásadě provádí n-krát operace plus resp. porovnání, které by měli být stejně
rychlé).
Úkol: Pomocí benchmarku ilustrujte, že třídění nemá lineární časovou složitost. Své pozorování ověřte prostřednictvím dokumentace k funkci
numpy.sort.
Funkce numpy.sort používá algoritmus quicksort, který má typickou časovou složitost
$O(n\cdot \log(n))$ (tu lze při použití pole s náhodným obsahem předpokládat s téměř
absolutní jistotou). Pro zobrazení této závislosti je nutné, aby se výrazně měnil
řád $n$.
from algbench import *
import numpy as np
r = benchmark(
np.sort, # při použití jedné funkce lze seznam vynechat
lambda n: np.random.uniform(0.0,1.0,n), (2**i for i in range(10,24)))
qplot(r, ["numpy.sort"])
2. Obecné kolekce¶
Obecné kolekce jsou datové struktury, které slouží k ukládání dat (tj. obecně objektů) a jsou využitelné v různých scénářích (tj. nikoliv jen pro určitý úzce definovaný účel). Obecné kolekce jsou často podporovány přímo běhovou podporou jazyka (aby byla jejich implementace rychlá a efektivní) resp. jsou naprogramovány prostřednictvím kódu v nízkoúrovňových jazycích (typicky je to jazyk C)
Obecné kolekce se mohou lišit:
- charakterem ukládaných hodnot
- rozhraním tj. nabídkou různých metod
- efektivitou jednotlivých operací (efektivita jednotlivých operací se může u jednotlivých kolekcí výrazně lišit)
Ve skutečnosti lze najít pouze jedinou operaci, která je sdílena všemi kolekcem, a tou je získání počtu položek (= uložených objektů).
Lineární kolekce¶
Lineární kolekce uchovávají uspořádanou posloupnost objektů, tj. lze je procházet v určitém neměnném pořadí. Implementují tak uspořádané množiny (ordered set) známé z matematiky.
Většina druhů lineárních kolekcí podporuje efektivní získání i-té položky tzv. indexaci (tento druh lineárních kolekcí se v Pythonu označuje jako sekvence).
Seznam¶
Základní lineární kolekcí jsou v Pythonu seznamy (list). Pokud seznamy neznáte, pak si si přečtěte stručný úvod do pythonských seznamů na http://diveintopython3.py.cz/native-datatypes.html (sekce 2.4).
Seznamy mají tuto charakteristiku:
- seznam obsahuje sekvenčně uložené odkazy na jednotlivé položky (nejsou uložené přímo v seznamu)
- indexace má konstantní časovou složitost
- do seznamu lze přidávat, ale efektivní je pouze přidávání na konec (časová složitost $O(1)$), přidávání na začátek nebo dovnitř seznamu má složitost $O(n)$
- prvky lze seznamu i odebírat (efektivní je opět jen odebírání na konci)
Python podporuje i neměnné sekvence, které se označují jako n-tice (tuple), ty jsou však užívány jen pro dočasné ukládání malého počtu prvků (maximálně malé desítky) a hrají tak v algoritmizaci jen pomocnou roli (nejčastěji se používají pro vracení více hodnot z funkcí).
Oboustranná fronta¶
Standardní knihovna Pythonu však podporuje i jednu dodatečnou sekvenci: oboustrannou frontu (deque), což je de facto seznam, u něhož který podporuje rychlé přidávání na obou koncích seznamu.
Většina implementací oboustranné fronty umožňuje indexaci s konstantní složitostí (i když často pomalejší než u seznamů).
Popis rozhraní pythonské oboustranné fronty viz https://docs.python.org/3/library/collections.html#collections.deque
Existuje několik možných implementací oboustranných front, z nichž rychle indexovatelné jsou implementace založené na po částech souvislých seznamech, detaily viz https://en.wikipedia.org/wiki/Double-ended_queue.
Úkol: pomocí benchmarkingu ověřte časovou charakteristiku oboustranné fronty (otestujte přidání na oba konce obou kolekcí)
import collections
from algbench import *
def list_append_right(n):
c = []
for i in range(n):
c.append(i)
def list_append_left(n):
c = []
for i in range(n):
c.insert(0, i)
def deque_append_right(n):
c = collections.deque()
for i in range(n):
c.append(i)
def deque_append_left(n):
c = collections.deque()
for i in range(n):
c.append(i)
r = benchmark([list_append_right, list_append_left, deque_append_right, deque_append_left],
init= lambda n: n, nrange=nrange(300))
qplot(r, ["list_right", "list_left", "deque_right", "deque_left"])
Výsledek je zřejmý. Přidávání na začátel seznamu (list_left)
je zřetelně nejpomalejší a má kvadratickou tendenci.
Oboustrannou frontu se při častém vkládání na začátek
vyplatí používat už i pro relativně malé seznamy (navíc ani vkládání
na konec je pro tuto velikost o trochu rychlejší).
Pole (arrays)¶
Další velmi důležitou sekvencí je pole. Polem se povětšinou rozumí lineární datová struktura, jež je optimalizována na rychlý přístup položkám prostřednictvím pozičního indexu, nepodporuje však přidávání na konec (pojem pole a seznam však obecně splývají a v některých jazycích mohou být odlišeny i jinou charakteristikou).
Python pole přímo nepodporuje jsou však implementovány v knihovně NumPy.
Pole v NumPy maji tyto charakteristiky:
1) podporují vícedimenziální indexaci (tj. pole representuje matematické objekty označované jako tenzory) 2) podporují přímé ukládaní hodnot číselných typů v položkách pole (seznam i v případě elementárních typů využívá odkazy na číselné objekty) 3) jsou kompatibilní s poli v nízkoúrovňových jazycích (C nebo Fortran). Nejsou sice identické (nesou běhové metainformace), ale datová část je organizována stejně.
Rozdíl mezi paměťovou implementací mezi pythonským seznamem a
Numpy polem ukazuje následující obrázek (převzatý z knihy Python Data Science Handbook).
Pole jsou klíčovou datovou strukturou pro datovou vědu a i my je budeme často využívat. Přehledný popis numpy polí najdete například na stránkách naucse,python.cz (https://naucse.python.cz/lessons/intro/numpy/). Prostudujte také originální dokumentaci na stránkách https://numpy.org.
Zaměřit byste se měli především na funkce vytvářející pole (včetně polí
vyplněná náhodnými hodnotami], indexaci včetně výřezů (slice), a skládání polí
pomocí funkcí vstack nebo hstack resp. pomocí r_ a c_ indexace
(je to tak trochu zneužití syntaxe indexace, ale je to užitečné).
Úkol:
Vytvořete následující pole (bez použití cyklů):
- jednorozměrné pole 100 prvků typu int8 naplněné nulami
- stejné pole ale naplněné čísly 0 až 99 (uložte je do proměnné
v)- stejné pole ale naplněné číslem 255
- stejné pole vyplněné náhodnými čísly v rozsahu 0 ... 255 s rovnoměrným rozdělením pravděpodobnosti
- jednorozměrné pole typu float naplněné čísly přesně padesáti ekvidistaními čísly od 100 do 200 (včetně)
- matici (= dvojrozměrné pole) položek typu double (float64) rozměru 5×5 bez vyplnění
- stejnou matici, ale vyplněné náhodnýmí čísly s normálním rozdělením se střední hodnotou 100. a rozptylem 100
- jednotkovou matici float64 čísel s rozměrem 100×100
- trojrozměrné pole 3×10×10 s jehož první vrstva je tvořena čísly 1.0, druhá 2.0 a třetí 3.0
- dvojrozměrné pole rozměrů 5×5, které je od levého horního rohu vyplněno po řádcích hodnotami 0.0 až 24.0 (uložte je do proměnné
m)
import numpy as np
v =np.arange(0,100, dtype=np.int8)
m = np.arange(0.0, 25.0).reshape(5,5)
print(v)
print(m)
Úkoly:
- z jednorozměrného vektoru
v(viz předchozí úkol) vytvořete výřez obsahující všechny prvky kromě prvního.- z vektoru
vvytvořete výřez obsahující pouze prvky se sudým indexem (počítáno od nuly)- z vektoru
vvytvořte výřez obsahující jen položky s lichými čísly (řešení musí fungovat pro pole s libovolným celočíselným obsahem)- z matice 'm' indexujte první dva řádky
- z matice 'm' indexujte poslední dva sloupce
- z matice indexujte podmatici tvořenou druhým a třetím řádkem a dtuhým až čtvrtým sloupcem
Otázky:
- V čem se liší paměťová representace výsledku úkolu 2 a 3?
- V čem se liší zápis
m[1:2, 2:3] od zápisum[1,2]`?
Spojový seznam¶
I přes bohatou nabídku kolekci v Pythonu v Pythonu a jeho standardních knihovnách, jedna z klasických lienárních kolekcí chybí: spojový seznam.
Spojový seznam (linked list) je tvořen jednoduchými datovými strukturami nazývanými uzly, které nesou hodnoty jednotlivých položek (tj. přímo jednotlivé položky nebo odkazy na ně) a vzájemně jsou spojeny sítí odkazů (refrencí). v případě jednosměrného seznamu je to vždy odkaz na následující uzel v v případě obosuměrného na předchozí.
Uzel, který nemá následovníka resp. předchůdce (u obousměrného seznamu) odkazuje
na objekt None. Zvnějšku je vždy dostupný první prvek resp. prvek poslední
(klíčové u obousměrného seznamu).
Mezi hlavní charakteristiky seznamu patří:
- konstantní čas přidávaní prvku kamkoliv do seznamu (je však nutno mít odkaz na prvek, za který resp. před který se prvek vkládá)
- neefektivní indexace (je nutno projít všechny prvky před hledaným, tj. časová složitost je průměrně $O(n)$).
- uzly nemusí ležet v paměti souvisle (což mimo jiné umožňuje flexibilnější úklid paměti)
Základní implementace jednosměrného seznamu je triviální, stačí nadefinovat třídu, jejiž instance jsou nové uzly.
%%writefile linked_list.py
class Node:
def __init__(self, val):
self.value = val
self.next = None # odkaz na počátku odkazuje na `None`
def insertNodeAfter(self, nextnode):
oldnext = self.next # uložíme si aktuální uzel
self.next = nextnode # atribut `next` nyní odkazuje na vkládaný uzel
nextnode.next = oldnext # a nastavíme odkaz na následující u vkládaného
return nextnode
def iter_to_end(self):
node = self
yield node.value
while node.next is not None:
node = node.next
yield node.value
Použití implementovaného spojového seznamuje zřejmé.
from linked_list import Node
start = Node(5)
start.insertNodeAfter(Node(6)).insertNodeAfter(Node(7))
start.insertNodeAfter(Node(5.5))
print(list(start.iter_to_end()))
Úkol: Rozhraní třídy
Nodeje bohužel trochu matoucí. MetodainsertNodeAftervkládá jediný uzel (pokud je jí předán spojový seznam, pak vloží jen jeho první prvek a zbytek vkládaného seznamu zahodí). To odpovídá běžnému chápání abstraktního rozhraní může to však být matoucí. Například pro postupné řetězení nelze využít následující kód.
chain = Node(5).insertNodeAfter(Node(6).insertNodeAfter(Node(7)))
print(list(chain.iter_to_end()))
Implementujte metodu
insertChainAfter, která umožní vkládat celé dílčí seznamy. Jaká je časová složitost této metody. Jakou časovou složitost má vytvoření seznamu o $n$ prvcích prostřednictvím vnitřního řetězení (tj. k prvnímu uzlu je připojen zbytek seznamu, který je opět vytvořen připojením zbytku seznamu k druhému uzlu, atd.
Časová složitost metody je $O(m)$, kde $m$ je délka přidávaného podseznamu (nezáleží na délce seznamu, do něhož se přidává, který je standardně označován jako $n$).
Časovou složitost vnitřního řetězení lze odvodit úvahou, že nejdříve se vytvoří poslední prvek, ten se přidá k předposlednímu (což si vyžádá traverzování přes jeden odkaz), výsledný dvouprvkový seznam se přidá k předpředposlednímu prvku (traverzování přes dva odkazy) a nakonec se seznam o dělce $n-1$ prvků přidá za první (s $n-1$ traverzováním). Celkový počet traverzovaní (každé vyžaduje porovnání a přiřazení) je tedy $1 + 2 + \cdots + (n-1) = \sum_{i=1}^{n-1} i$, což je součet jednoduché arimetické posloupnosti tj. $\frac{(n-1)\cdot n}{2}$. Časová složitost je tedy řádu $O(n^2)$, je tedy podstatně horší než vnější řetězení uvedené výše (to vkládá prvky od prvního do posledního) se složitosti $O(n)$. Cesta do pekel je dlážděna dobrými úmysly.
Datové struktury optimalizované pro vyhledávání¶
Datové struktury optimalizované pro vyhledávání typicky ukládají dvojice objektů označovaných jako klíč a hodnota, přičemž nabízejí rychlou implementaci operace vyhledávání podle klíče (s asymtotickou časovou složitostí lepší než lineární). Rychlé by měly být i operace vkládání i výmazu (především pokud se předpokládá, že bude docházet k časté aktualizaci datové struktury).
Hashovací tabulky¶
Hashovací tabulky jsou nejčastěji používanými vyhledávacími datovými strukturami, které typicky dosahují ve všech třech základních operacích (vložení, vyhledávání, výmaz typické asymptotické složitosti $O(1)$, tj. čas prakticky nezávisí na velikosti kolekce.
Toto chování je však dosaženo jen při splnění následujících podmínek:
- existuje hashovací funkce, která libovolný vložený klíč zobrazí na celé číslo v rozsahu $0\cdots m$, tak že zajistí prosté zobrazení do množiny (tj. každý vstupní klíč zobrazí na jiné číslo) nebo minimalizuje počet kolizí (situací, kdy jsou dva či více klíčů zobrazeny na stejnou hodnotu). Průměrný počet kolidujících klíčů by neměl růst s velikostí kolekce (tj. počtem klíčů).
je k dispozici dostatečné množství paměti, které závisí primárně na hodnotě $m$, ale i na $n$ (pokud je $n > m$)
Hashovací tabulky můžete prostudovat na Wikipedii (článek na anglické Wikipedii je instruktivní a obsahuje přehled existujích přístupů), resp. (přímo pro Python, pak doporučuji pythonds.
Pozornost věnujte především problematice řešení kolizí (řetězení, otevřené násobení) a důvodům pro tzv. prehashování.
Hashovací tabulky jsou v Pythonu použity pro implementaci slovníků (dict)
a množin (set, frozenset).
Implementace slovníků využívá otevřené adresování, kde se jako sondovací sekvence (probing sequence) používá jednoduchý generátor náhodných čísel. Velikost tabulky je mocninou dvou (počínaje osmi položkami) a tabulka je zdvojnásobena a rehashována pokud její zaplnění vzroste nad 2/3. Nověji se hashování provádí nad tabulkou indexů, které odkazují do lineární tabulky (viz příspěvek na stackoverflow).
Otázka: Proč je novější implementace slovníků (od verze 3.6) používá nepřímé adresování prostřednictvím seznamu indexů? Diskutujte především použití slovníků pro předávání pojmenovaných parametrů.
Python podporuje obecně i podporu hashovacích funkcí. Téměř všechny standardní
třídy u nichž to dává smysl (instance jsou neměnné a lze testovat jejich
shodu) jsou v Pythonu hashovatelné tj. implementují speciální metodu __hash__.
Tato obecná metoda vrací číslo v rozmezí 0 až maxint. Reálně použitá
hashovací funkce volá obecnou hashovací metodu, využívá však jen příslušný
počet bitů (nižších) bytů.
Tyto obecné hashovací metody, mají relativně dobré chování pro náhodné vstupní hodnoty resp. hodnoty z běžně používaných podmnožin (např. hashování řetězcových hodnot má dobré chování i pro texty ve většině světových jazyků v běžných kódováních). Pro specializované skupiny řetězců však může být podíl kolozí větší.
V případě uživatelských tříd, u nichž se předpokládá, že busdou sloužit jako klíče, resp. budou ukládány do množiny je nutno definovat vlastní speciální metodu hash.
Tato metoda by měla splňovat tyto podmínky:
- výsledek by měl refletovat všechny dílčí hodnoty a jejich bity
- i malá změna dat (např. záměna jednoho písmene) by měla vést ke změně výsledné hodnoty
- hashované hodnoty by měli využívat všechny bity výsledku (výsledek je 64 bitový), pokud to není možné tak alespoň $n$-bitů, kde je $n$ dostatečně velké
- výpočet by měl být rychlý nejlépe s konstantní dobou provedení (nezávislou na rozsahu vtsupnách hodnot)
Úkol: Diskutujte následující standardní implementaci metody
__hash__(zkráceno a převedeno do Pythonu, originál je v jazyce C), vzhledem k výše uvedeným požadavkům.value = ord(self[0]) << 7 for char in self: value = (1000003 * value) ^ ord(char) # násobení je omezeno na 64 bitů value = value ^ len(self) if value == -1: value = -2 return value
Pro složené objekty se preferuje využití operace xor (zápis v Pythonu ^)
mezi všemi datovými členy, které se využívají v operaci testování shody.
To je ve většíně případů dobrá volba (operace xor je rychlá zachovává
rovoměrné rozložení bitů, nemusí však být vždy optimální,
viz například následující úkol):
Úkol Implementujte třídu pro representaci pozice šachové figurky využívající atributy
rowacolumn. Implementujte speciální metodu pro operátor rovnosti (__eq__) a obecnou hashovací metodu. Diskutujte metodu__hash__a pokuste se najít nejlepší řešení.
class ChessPosition:
def __init__(self, row, column):
self.row = row
self.column = column
def __eq__(self, other):
return self.row == other.row and self.column == other.column
def __hash__(self):
return hash(self.row) ^ hash(self.column)
p1 = ChessPosition(1,1)
p2 = ChessPosition(2,2)
print(hash(p1))
print(hash(p2))
Hashovací funkce využívající standardní přístup vrací pro všechny pole na diagonále
(operace xor aplikovaná na dvě shodné hodnoty vrací vždy nulu). Jinak řečeno 1/8
možných instancé této třídy vrací stejnou hashovací hodnotu.
O něco lepším řešením by byla hashovací hodnota row*8 + column, která však v nižších
bytech zohledňuje jen sloupec a nikoliv řádek (tj. při velikosti hashovací tabulky například
8 bytů rozhoduje spadají všechny pole nav jedmom sloupci do stejného místa tabulky,
což však není překvapivé).
Binární vyhledávací stromy¶
Binární vyhledávací stromy (binary search tree) jsou klasickou datovou strukturou, která je založena na jednoduchém a elegantním přístupu. V typickém případě mají základní oprace nad binárními stromy časovou složitost $O(log(n))$, což je sice horší než v případe hashovacích tabulek, pro běžné velikosti kolekcí (= kolekcí, které se vejdou do paměti) však nejsou rozdíly velké.
Binární vyhledávací stromy lze aplikovat pouze pro položky, pro něž je definováno
rozumné uspořádání tj. kromě relace rovnosti a nerovnosti
musí být definovány i operace <, <=, > a >=, tak aby mezi nimi platily
základní vztahy a uspořádání definované relací <= byli tzv. úplné (viz
https://en.wikipedia.org/wiki/Total_order).
Kompletní a dobře zdokumentovanou implementaci binárních stromů najdete na stránkách runestone.academy (https://runestone.academy/runestone/books/published/pythonds/Trees/BinarySearchTrees.html).
Následující kód je převzat z výše uvedené implementace, je mírně upraven a okomentován:
%%writefile binarysearchtree.py
import collections.abc
class TreeNode: # pomocná třídy: uzel binárního stromu
def __init__(self,key,val,left=None,right=None,parent=None):
self.key = key
self.payload = val
self.leftChild = left
self.rightChild = right
self.parent = parent
def hasLeftChild(self):
return self.leftChild is not None
def hasRightChild(self):
return self.rightChild is not None
def isLeftChild(self):
return self.parent and self.parent.leftChild == self
def isRightChild(self):
return self.parent and self.parent.rightChild == self
def isRoot(self):
return not self.parent
def isLeaf(self):
return not (self.rightChild or self.leftChild)
def hasAnyChildren(self):
return self.rightChild or self.leftChild
def hasBothChildren(self):
return self.rightChild and self.leftChild
def spliceOut(self): # pomocná metoda pro výmaz
if self.isLeaf():
if self.isLeftChild():
self.parent.leftChild = None
else:
self.parent.rightChild = None
elif self.hasAnyChildren():
if self.hasLeftChild():
if self.isLeftChild():
self.parent.leftChild = self.leftChild
else:
self.parent.rightChild = self.leftChild
self.leftChild.parent = self.parent
else:
if self.isLeftChild():
self.parent.leftChild = self.rightChild
else:
self.parent.rightChild = self.rightChild
self.rightChild.parent = self.parent
def findSuccessor(self): # pomocná metoda pro výmaz
succ = None
if self.hasRightChild():
succ = self.rightChild.findMin()
else:
if self.parent:
if self.isLeftChild():
succ = self.parent
else:
self.parent.rightChild = None
succ = self.parent.findSuccessor()
self.parent.rightChild = self
return succ
def findMin(self): # pomocná metoda pro výmaz
current = self
while current.hasLeftChild():
current = current.leftChild
return current
def replaceNodeData(self,key,value,lc,rc): # # pomocná metoda pro výmaz
self.key = key
self.payload = value
self.leftChild = lc
self.rightChild = rc
if self.hasLeftChild():
self.leftChild.parent = self
if self.hasRightChild():
self.rightChild.parent = self
def __iter__(self): # iterátor přes (seřazené) klíče podstromu
if self:
if self.hasLeftChild():
for elem in self.leftChild:
yield elem
yield self.key
if self.hasRightChild():
for elem in self.rightChild:
yield elem
class BinarySearchTree(collections.abc.MutableMapping):
def __init__(self):
self.root = None
self.size = 0
def length(self):
return self.size
def __len__(self):
return self.size
def __iter__(self):
return self.root.__iter__()
def put(self,key,val):
if self.root:
self._put(key,val,self.root)
else:
self.root = TreeNode(key,val)
self.size = self.size + 1
def _put(self,key,val,currentNode):
if key < currentNode.key:
if currentNode.hasLeftChild():
self._put(key,val,currentNode.leftChild)
else:
currentNode.leftChild = TreeNode(key,val,parent=currentNode)
else:
if currentNode.hasRightChild():
self._put(key,val,currentNode.rightChild)
else:
currentNode.rightChild = TreeNode(key,val,parent=currentNode)
def __setitem__(self,k,v):
self.put(k,v)
def get(self, key, default):
if self.root:
res = self._get(key,self.root)
if res:
return res.payload
else:
return default
else:
return default
def _get(self,key,currentNode):
if not currentNode:
return None
elif currentNode.key == key:
return currentNode
elif key < currentNode.key:
return self._get(key,currentNode.leftChild)
else:
return self._get(key,currentNode.rightChild)
def __getitem__(self,key):
if self.root:
res = self._get(key,self.root)
if res:
return res.payload
else:
raise KeyError("Invalid key")
else:
raise KeyError("Invalid key")
def __contains__(self,key):
if self._get(key,self.root):
return True
else:
return False
def delete(self,key):
if self.size > 1:
nodeToRemove = self._get(key,self.root)
if nodeToRemove:
self.remove(nodeToRemove)
self.size = self.size-1
else:
raise KeyError('Error, key not in tree')
elif self.size == 1 and self.root.key == key:
self.root = None
self.size = self.size - 1
else:
raise KeyError('Error, key not in tree')
def __delitem__(self,key):
self.delete(key)
def remove(self,currentNode):
if currentNode.isLeaf(): #leaf
if currentNode == currentNode.parent.leftChild:
currentNode.parent.leftChild = None
else:
currentNode.parent.rightChild = None
elif currentNode.hasBothChildren(): #interior
succ = currentNode.findSuccessor()
succ.spliceOut()
currentNode.key = succ.key
currentNode.payload = succ.payload
else: # this node has one child
if currentNode.hasLeftChild():
if currentNode.isLeftChild():
currentNode.leftChild.parent = currentNode.parent
currentNode.parent.leftChild = currentNode.leftChild
elif currentNode.isRightChild():
currentNode.leftChild.parent = currentNode.parent
currentNode.parent.rightChild = currentNode.leftChild
else:
currentNode.replaceNodeData(currentNode.leftChild.key,
currentNode.leftChild.payload,
currentNode.leftChild.leftChild,
currentNode.leftChild.rightChild)
else:
if currentNode.isLeftChild():
currentNode.rightChild.parent = currentNode.parent
currentNode.parent.leftChild = currentNode.rightChild
elif currentNode.isRightChild():
currentNode.rightChild.parent = currentNode.parent
currentNode.parent.rightChild = currentNode.rightChild
else:
currentNode.replaceNodeData(currentNode.rightChild.key,
currentNode.rightChild.payload,
currentNode.rightChild.leftChild,
currentNode.rightChild.rightChild)
Výše uvedená implementace binárního stromu implementuje abstraktní datový
typ collections.abc.MutableMapping
(viz https://docs.python.org/3/library/collections.abc.html) a to
prostřednictvím dědění (kromě přímo definovaných metod tak podporuje
i zděděné tzv. mixin metody jako je metoda keys, items, values, i když
ne příliš efektivně).
from binarysearchtree import BinarySearchTree
tree = BinarySearchTree()
# vložení klíčů a hodnot
tree[2] = "dve"
tree[5] = "pet"
tree[4] = "ctyri"
tree[1] = "jedna"
tree[3] = "tri"
# nalezení hodnoty podle klíče
print(tree[2])
# iterace
for key in tree.keys():
print(key)
Úkol: Nakreslete tvar výsledného binárního stromu (stačí uvést pouze klíče).
Rada: Zkuste použít nějaký vizualizátor binárních vyhledávacích stromů (např. http://btv.melezinek.cz/binary-search-tree.html)
Na rozdíl od minimalistické implementace spojového seznamu, implementuje výše uvedený kód obalovou třídu (tj. nikoliv jen stavební uzel). Kód je přehledný a samopopisný. Výjimkou je kód pro výmaz položky, který je dlouhý a navíc používá netriviální pomocné metody. Zde je detailnější prostudování implementace nezbytné.
Shrňme si základní charakteristiky jednoduchých binárních vyhledávacích stromů:
- typická časová složitost všech základních operací je $O(log(n))$, operace výmazu je však díky komplexnější implementaci méně efektivní (především u menších kolekcí).
- přes kolekci lze iterovat tak, že pro každý prvek platí, že je větší nebo rovný něž prvek předcházející (samozřejmě kromě prvního), tj. posloupnost je seřazená.
Triviální binární vyhledávací stromy však mají i problematické vlastnosti.
- rychlost všech operací závisí na hloubce stromu (= maximální vzdálenost od kořene k listu), pro kterou platí $h\geq \lfloor \log _{2}n\rfloor$. Pokud je stromu hloubka rovna minimální hodnotě (tj. $\lfloor \log _{2}n\rfloor$) pak se strom označuje jako vyvážený. Vyvážené stromy optimalizují rychlost všech základních operací.
V praxi dostačují i tzv. téměř vyvážené stromy, jejichž hloubka je rovna $k\cdot\lfloor \log _{2}n\rfloor$, kde $k \ll n$ (tj. je to dostatečně malé číslo, typicky menší než 3). Operace v téměř vyvážených stromech jsou pomalejší jejich asympotická časová složitost je však stále rovna $O(log(n))$.
Pokud jsou položk vkládány v náhodném pořadí, je strom téměř vždy téměř vyváženy. Pokud však čísla vkládáme již seřazená, pak vznikne absolutně nevyvážený strom s jednou dlouhou větví (bez větvení). Tento strom má hloubku $n$ (má jen jeden list) a jedná se de facto o spojový seznam, v němž mají všechny základní operace časovou složitost $O(n)$.
Úkol: Implementujte funkci, která pro daný strom zjistí jeho hloubku. Ověřte, že hloubka stromu vzniklého z náhodné posloupnosti klíčů je téměř optimální.
from binarysearchtree import BinarySearchTree, TreeNode
from numpy.random import rand
def tree_depth(tree: BinarySearchTree):
return _depth(tree.root)
def _depth(node: TreeNode):
if node is None:
return 0
return 1 + max(_depth(node.leftChild), _depth(node.rightChild))
def random_tree(n):
tree = BinarySearchTree()
for i,x in enumerate(rand(n)):
tree[x] = i
return tree
for _ in range(5):
t = random_tree(2**20-1) # cca 1M prvků (optimální hloubka 20)
print(tree_depth(t))
Empirické výsledky ukazují, že se hloubka stromu pohybuje kolem 50, což je sice více než hloubka optimální (ta je 20), je to však stále o několik řádů méně než počet počet prvků.
Úkol: Ověřte rychlost vytváření a vyhledávaní v binárních vyhledávacích stromech ve srovnání s hashovacími tabulkami pro různé řády hodnot n.
Samovyvažující se binární vyhledávací stromy¶
V případech, kdy nelze zaručit dostatečnou náhodnost vstupu (tj. například vstupní posloupnost obsahuje uspořádané podsekce), je nutné použít tzv. samovyvažujících se stromů. U těchto stromů dochází při vkládání k tzv. vyvažování tj. přeorganizování jednotlivých větví tak, aby byl strom téměř vyvážen (cílem není zcela optimální vyvážení).
Červenočerné stromy¶
Červenočerné stromy patří mezi nejoblíbenější samovyvažující se stromy, nabízejí rozumné vyvážení (maximální hloubka je dvojnásobek optimální) při relativně pochopitelné sémantice a algoritmech.
Detailní popis červenočerných stromů naleznete na stránkách https://algorithmtutor.com/Data-Structures/Tree/Red-Black-Trees/ a to včetně implementace.
Úkol: Porovnejte výkonnost červenočerného a běžného binárního vyhledávacího stromu na náhodných datech resp. datech uspořádaných.
AVL stromy¶
Dalšími běžně používanými samovyvažujícími se stromy jsou AVL stromy. Detailní popis včetně implementace najdete na https://runestone.academy/runestone/books/published/pythonds/Trees/AVLTreeImplementation.html.
AVL stromy jsou de facto speciálním případem červenočerných stromů, dosahují však teoreticky lepšího vyvážení. Jejich hloubka je v nejhorším případě 0,72 krát lepší než je tomu u nejhoře vyvážených červenočerných stromů). Vyžadují však více uložených dat a implementace základních operací je o něco složitější.
Úkol: Porovnejte hloubku AVL a
3. Speciální kolekce¶
Speciální kolekce se v zásadě neliší od dříve uvedených kolekcí univerzálních, většinou však mají specifičtější rozhraní a neslouží tudíž primárně jako úložiště dat.
Orientované grafy¶
Python ve své standardní podobě nepodporuje representaci orientovaných grafů resp. odvozených struktur (neorientovaných grafů, multigrafů, ohodnocených grafů apod.)
Jednotlivé typy grafů lze representovat mnoha způsoby například:
- maticí sousednosti
- řídkou maticí sousednosti
- seznamem hran
- slovníkem, jehož klíčem je označení uzlu a hodnotou seznam
- slovníkem, jehož klíčem je označení uzlu a hodnotou seznam sousedních vrcholů
Na stránce https://www.python-course.eu/graphs_python.php najdete detailní popis implementace třídy pro neorientované grafy.
Otázka: Výše uvedená implementace je prezentována jako implementace neorientovaných grafů, ve skutečnosti je však vhodnější pro representace grafů orinetovaných. Uveďte proč?
Pro praktickou práci s grafy je však vhodnější využít komplexnějších knihovnu. Dporučuji knihovnu NerworkX. Prostudujte tutoriál knihovny: NetworkX tutorial a také první kapitolu knihy NetworkX Reference resp. podle vašich zájmů i další části referenční příručky.
Otázka: Knihovna používá pro representaci slovník slovníků (slovníků). Jaké výhody tato implementace přináší (paměťová a časová efektivita) a jaké jsou její nevýhody (omezení).
Úkol: Vytvořte úplný neorientovaný graf $K_5$, jehož uzly jsou číslovány 0 až 5. Zobrazte ho. Zkuste nízkoúrovňový přístup i vestavěný generátor.
!pip3 install networkx
Řešení pomocí nízkoúrovňového přídávání hran (zdrojem je iterátor přes všechny dvojčlenné kombinace).
import networkx as nx
import itertools
import matplotlib.pyplot as plt
g = nx.Graph()
for i in range(10):
g.add_node(i)
g.add_edges_from(itertools.combinations(range(10),2))
nx.draw(g)
plt.show()
Řešení pomocí vestavěného generátoru.
g = nx.complete_graph(10)
nx.draw_shell(g)
plt.show()
Úkol: Vytvořte neorientovaný ohodnocený graf, representující níže uvedenou tabulku dojezdových minimálních časů mezi krajskými městy (územních) krajů. Graf zobrazte.
BR CB HK PR BR CB 02:21:00 HK 02:14:00 OS 01:42:00 PR 01:58:00 01:40:00 01:16:00 PL 02:00:00 01:02:00 UL 00:59:00
d = nx.Graph()
d.add_edge("CB","BR", weight="2:21")
d.add_edge("HK","BR", weight="2:14")
d.add_edge("OS","BR", weight="1:42")
d.add_edge("PR","BR", weight="1:58")
d.add_edge("PR","CB", weight="1:40")
d.add_edge("PR","HK", weight="1:16")
d.add_edge("PL","CB", weight="2:00")
d.add_edge("PL","PR", weight="1:02")
d.add_edge("UL","PR", weight="0:59")
for f,t,w in d.edges(data="weight"):
h,m = w.split(":")
d[f][t]["weight"] = int(h) * 60 + int(m)
nx.draw(d, with_labels=True)
plt.show()
Vykreslování grafů je překvapivě komplexní problematika, kde nasazení běžných obecných algoritmů nevede k příliš kvalitním výsledkům. NetworkX využívá imlicitně tzv.pružinový algoritmus (viz Force-directed graph drawing), jehož výsledky závisí na počáteční (náhodné) inicializaci. Je proto potřeba využít opakované vyhodnocení příslušného kódu. Stabilnější výsledky poskytuje algoritmus umisťující uzly do kruhů či soustředných kruhů (draw_circular, draw_shell), při tomto přístupu, však mohou vznikat i zbytečná křížení hran. Použít lze i zcela náhodné umisťování uzlů (draw_random).
Fronty a zásobníky¶
Fronty a zásobníky jsou sekvenční kolekce s rozhraním, které je v zásadě omezeno jen na dvě operace: vkládání položky a její vyjímání. Běžně bývá podporováno i testování prázdnosti. Obě speciální kolekce se bežně implementují jako adaptéry nad obecnými sekvencemi (tj. interně využívají obecnější sekvence, navenek však poskytují jen obě výše uvedné metody resp. vlastnost).
Oba klasické abstraktní datové typy se liší sémantikou operace vyjímání. V případě zásobníků je vyjmuta poslední vložená (a nevyjmutá) položka v případě fronty je vyjmuta položka nejdřívě vložená (a samozřejmě ještě nevyjmutá).
Základní anglickou terminologii spojenou s frontami a zásobníky ukazuje následující tabulka:
| fronta | zásobník | |
|---|---|---|
| zákl. označení | queue [kjuː] | stack [stæk] |
| mechanismus | First In First Out | Last In Last Out |
| zkratka | FIFO | LIFO |
| operace In | enqueue | push |
| operace Out | dequeue | pop |
Zásobníky jsou v Pythonu snadno implementovatelné pomocí běžného seznamu a operací append a pop (s časovou složitostí $O(1)$). Pro testování prázdnosti lze využít běžný boolovský test (seznam je "pravdivý" tehdy pokud je neprázdný).
Úkol: Naimplementujte adaptér nad seznamem. který bude podporovat jen základní operace zásobníku (testování prázdnosti pomocí speciální metody
__bool__).
class Stack:
def __init__(self):
self.data = []
def __bool__(self):
return bool(self.data)
def push(self, item):
self.data.append(item)
def pop(self):
return self.data.pop()
# test
s = Stack()
for i in range(5):
s.push(i)
while(s): # dokud je neprázdný
print(s.pop())
Čísla se vypisují v opačném gardu (první je vyjmuto a vypsáno poslední vložené).
Frontu lze teoreticky také implementovat pomocí seznamu, jen je nutné prvky vkládat na začátek seznamu. To je však neefektivní, neboť časová složitost této implementace by byla $O(n)$, kde n v tomto případě označuje aktuální délku fronty (tj. pokud je průměrná či ještě l=pe maximální aktuální délka fronty relativně malá v řádu malých desítek, pak je tato implementace bez problémů použitelná).
Stejnou efektivitu s konstatní časovou složitostí však nabízí obostranná fronta, u níž se využije jen vkládání na začátek a vyjímání na konci.
Zásobníky jsou v Pythonu snadno implementovatelné pomocí běžného seznamu a operací append a pop (s časovou složitostí $O(1)$). Pro testování prázdnosti lze využít běžný boolovský test (seznam je "pravdivý" tehdy pokud je neprázdný).
Úkol: Naimplementujte adaptér nad oboustrannou frontou, který bude podporovat jen základní operace fronty (a testování prázdnosti pomocí speciální metody
__bool__).
import collections
class Queue:
def __init__(self):
self.data = collections.deque()
def __bool__(self):
return bool(self.data)
def enqueue(self, item):
self.data.appendleft(item)
def dequeue(self):
return self.data.pop()
# test
s = Queue()
for i in range(5):
s.enqueue(i)
while(s): # dokud je neprázdný
print(s.dequeue())
Pythonská standardní knihovna obsahuje i modul queue, která podporuje synchronizované fronty použitelné pro vícevláknové programy (pro programy víceprocesorové existuje implementace synchronizačních front v modulu
multiprocessing).
Úkol: Prostudujte dokumentaci na adrese https://docs.python.org/3/library/queue.html.
Synchronizační fronty se od běžných front liší v následujících charakteristikách:
- lze je využít pro mezivláknovou či meziprocesorovou komunikaci
- operace
get'odpovídající operacienqueuemůže být tzv. blokující tj. vlákno/proces čtoucí z fronty je pozastaven, dokud jiné vlákno-proces do fronty nevloží další položku. - fronta může mít nastavenu maximální velikost, pokud je nastavena, tak se může blokovat i vkládací operace
put. - fronty lze využít i v jednovláknovém programu, nikdy však nesmí dojít k blokujícímu volání (tj. např. nelze vyjímat prvek z prázdné fronty).
Kromě klasické (a nejužitečnější) fronty je implementována i tzv. LifoQueue (trochu
podivné kontradiktorické jméno je zde využito i pro zásobník) a tzv. prioritní fronta
PriorityQueue, což je speciální kolekce se sémantikou podobná frontě, ale s možností předbíhání
(prioritní fronta není frontou podle běžné definice fronty!).
Poznámka: Použití netradičních jmen jako je
LifoQueuea alternativních jmen je v některých kontextech dosti běžné. Při použití je proto důležité znát sémantiku každé kolekce a jejích metod.Úkol: Využijte priorotní frontu pro representaci fronty budoucích událostí v kalendáři. Každá událost je dána počátečním časem a prioritou 0 … 5 (5 je nejvyšší). Jako první je vyjímána událost, která nastane nejdříve, pokud je více takovýchto událostí, pak ta s vyšší prioritou . Rada: pro ukládání využijte trojice, v níž posledním členem je popis události resp. vlastní třídu s vhodně definovaným uspořádáním prvků.
Cyklické fronty (buffery)¶
Termínem cyklické fronty (circullar buffer) se označuje specifická implementace front. Cyklické fronty mají omezenou velikost, lze je však implementovat prostřednictvím polí (tj. souvislých sekvencí fixních velikostí)
Úkol: Nastudujte cyklické fronty (viz například hezký popis a animací na Wikipedii) a implementujte je jako třídu v Pythonu.
Cyklické fronty poskytují rychlé operace, jejiž časová složitost je i v nejhorším případě konstantní (není zde žádná dodatečná režie pro realokaci, apod.), oproti oboustranným frontám však mohou vyžadovat více paměti (alokována je vždy maximální velikost bez ohledu na uaktuální počet prvků ve frontě).
Řídké matice¶
Především pro strojové učení jsou klíčovými datovými strukturami tzv. struktury efektivně representující tzv. řídké matice či obecněji pole (sparse matrix(array)). Jako řídká matice (či obecněji řídká pole) se označují matice, u nichž napoloviční počet položek obsahuje nulovou hodnotu (obecněji jakoukoliv fixní hodnotu, i když tento typ řídkých matic se v praxi příliš nevyskytuje).
Řídké matice lze samozřejmě representovat pomocí běžných polí, je to však neefektivní.
import numpy as np
from random import randrange
def random_sparse_matrix(shape, values):
n = len(values)
m = np.zeros(shape, dtype=values.dtype)
while(n > 0):
x,y = randrange(shape[0]), randrange(shape[1])
if m[x,y] == 0:
m[x,y] = values[n-1]
n -= 1
return m
print(random_sparse_matrix((10,10), np.arange(0,11)))
Hustota této řídké matice je jen 10% (= obsahuje jen 10% nenulových prvků) jenak řečeno řídkost (sparsity) je 90%. Matice zaujímá v paměti 10×10x8 = 800 bytů v paměti (pro 64bitové systémy).
Jednotlivé representace řídkých matic jsou uvedeny na Wikipedii. Detailnější popis a především dokumentaci rozhraní naleznete na Scipy.
Otázka: Jaké z podporovaných formátů jsou vhodné pro konstrukci matice a jaké pro klíčové operace jako je násobení matic či výpočet inverzní matice?
Úkol: V praxi vyřešte soustavu rovnic zadanou pomocí řídké (avšak regulární) matice (bez konverze).
4. Algoritmy rozděl a panuj¶
Algoritmy typu rozděl a panuj (divide-and-conquer) tvoří důležitou třídu algoritmů, které pracují na principu rekurzivního dělení problémů na dva či více (jednodušších) podproblémů. Tyto podproblémy jsou v rámcu rekurze děleny na podproblémy druhé úrovně a to pokračuje tak dlouho dokud nezískáme podproblémy, které jsou dostatečně malé aby jejich řešení bylo triviální (např. u třídících algoritmů je triviálním problémem seřazení jednoprvkových nebo dvojprvkových sekvencí).
Přirozeným programátorským přístupem pro algoritmy tohoto druhu je rekurze. Rekurze však přináší dodatečné paměťové nároky a nemusí být dostatečně efektivní. Pro určitou třídu algoritmů ntohoto druhu však lze nalézt nerekurzivní řešení a to buď přímým přepisem (je u algoritmů s tzv. koncovou rekurzí, které však nejsou typickými algoritmy rozděl a panuj, neboť v každém kroku vznikne jen jedna podúloha, viz níže) resp. komplexnějším přístupem.
Algoritmy sniž-a-panuj¶
Tyto algoritmy jsou založeny nalezení jediné podúlohy, která výrazně sníží její rozsah (typické je například půlení, kdy je podúloha má nejvýše poloviční rozsah). Klasickou úlohou je binární vyhledávání, které bylo známo již ve starověku.
Binární vyhledávání¶
Binární vyhledávání je jednoduchý a přitom dostatečně efektivní algoritmus pro vyhledávání položky v seřazeném (setříděném) poli s časovou složitostí $O(n)$.
Popis i analýzu algoritmu najdete na https://runestone.academy/runestone/books/published/pythonds/SortSearch/TheBinarySearch.html.
Otázka: Jaké jsou výhody a nevýhody tohoto algoritmu ve srovnání s využitím specializovaných vyhledávacích kolekcí jako je binární strom nebo hashovací tabulka?
Binírní vyhledávání je základem standardního modulu bisect, který kromě vyhledávání podporuje i vkládání (hledá i pozice vhodné pro vložení), dokumentace viz https://docs.python.org/3.8/library/bisect.html.
Úkol: Implementujte binární pomocí rekurze i bez rekurze. Určete paměťovou složitost obou řešení a porovnejte jejich rychlost.
Efektivní řadící algoritmy¶
Nejefektivnější univerzální řadící algoritmy jsou založeny na mechanismu rozděl a panuj. Jedná se především o algoritmy quicksort a mergesort.
Popis algoritmu quicksort včetně jeho implementace najdete na https://runestone.academy/runestone/books/published/pythonds/SortSearch/TheQuickSort.html. Tento článek nediskutuje volbu tzv. pivota, což může výrazně ovlivnit efektivitu pro některé vstupy. Dobrým východiskem této klíčové diskuse je článek na wikipedii.
Úkol: Diskutujte výhody a nevýhody volby pivota mediánem (ze tří) a náhodné volby mediánu (nejhorší časová složitost, reálná rychlost).
Algoritmus mergesort, je ve většině případů pomalejší než quicksort, ale je stabilní (tj. zachovává pořadí položek se stejným řadícím klíčem) a i v nejhorším případě má časovou složitost $O(n log(n)$
Implementace algoritmu mergesort v Pythonu (včetně diskuse) je uvedena např. v
https://runestone.academy/runestone/books/published/pythonds/SortSearch/TheMergeSort.html (prostudujte).
Úkol: Porovnejte rychlost obou třídících algoritmů (u
quicksorts triviální volbou pivota a 3-mediánem ) pro pole s náhodným čísleným obsahem resp. již setříděná pole (včetně setříděných o opačné gardu). Interpretujte výsledky.
5. Hladové algoritmy¶
Jako hladové algoritmy se označují algoritmy, které hledají určité (globální) optimum tím, že provádějí lokální rozhodnutí na lokální úrovni.
Hladové algoritmy obecně nemusí najít optimální řešení a ve specifických případech mohou dokonce najít to nejhorší možné. Pro širokou třídu problémů však vedou k řešení, které je akceptovatelně blízké skutečnému a je tak využít jako heuristiky. Pro určité třídy problémů hladové algoritmy najdou vždy nejlepší řešení.
Třídě hladových algoritmů jsou blízké i algoritmy, jejichž cílem není optimalizace, ale prohledání určitého stavového prostoru či seřazení prvků založené na principu lokální volby.
Elementárním a v praxi často používaným algoritmem je algoritmus hledající nejmenší počet mincí či bankovek potřebných pro zaplacení určité sumy. V tomto algoritmu stačí v každém kroku zvolit vždy největší minci/bankovku menší než je zbývající suma (lokální rozhodnutí). Tento algoritmus nemusí obecně vést k optimálnímu rozhodnutí, v případě měnového systému v ČR (a mnoha dalších zemí) však vždy nalezne optimální řešení.
Obecné řešení vyžaduje využití dynamického programování viz například https://www.geeksforgeeks.org/find-minimum-number-of-coins-that-make-a-change/ nebo článek na Wikipedii.
Poznámka: Zjištění zda lze pro určité nominální hodnoty mincí bankovek najít optimální řešení pomocí hladového algoritmu je (tj. zda jde o tzv kanonický systém) není triviální, viz Xuan Cai and Yiyuan Zheng, nonical Coin Systems for Change-Making problems
Algoritmy poskytující optimální řešení¶
Grafové algoritmy¶
Důležitou skupinu hladových algoritmů tohoto typu jsou algoritmy pracující nad ohodnocenými grafy.
Klíčovou úlohou v ohodnocenných grafech hrají algoritmuy hledání nejkratší cesty, kde délka cesty je sumou ohodnocení všech hran na dané cestě. Tento algoritmy má široké nasazení v různých oblastech např, pro hledání optimálního routování apod.
Základním algoritmem v této oblasti je Dijkstrův algoritmus: viz například https://math.mit.edu/~rothvoss/18.304.3PM/Presentations/1-Melissa.pdf.
Úkol: Implementujte základní Dijkstrův algoritmus v pro NetworkX grafy a porovnejte se implementací v knihovně-
Alternativně lze využít v mnoha případech efektivnější algoritmus A* (viz Wikipedia a zajímavou presentaci na http://voho.eu/wiki/algoritmus-a-star/).
Úkol: Využijte vestavěnou implementaci algoritmu
A*v knihovně NetworkX. Použijte tento algoritmus pro nalezení nejkratší cesty (mezi dvěma vzdálenými) v síti, kterou budou tvořit vybrané geografické body spojené hranami, pokud jsou n-tými nejbližšími sousedy (měřeno přímou vzdálenosti), a hrany jsou ohodnoceny přímou vzdáleností. Jako heuristickou funkci zvolte přímou vzdálenost k cílovému bodu. Příklad: polohy obcí v ČR (jsou k dohledání na internetu viz např. https://github.com/33bcdd/souradnice-mest), okolí tvoří 4 nejbližší obce, najděte nejkratší cestu například mezi obcemi Rajhradice a Dobříň. Rada: uzly musí nést informace o geografických souřadnicích (např. atributylongitudealatitude), přímá vzdálenost mezi body by měla být representována atributemweight.
Zajímavou třídou hladových algoritmů, jsou algoritmy hledání tzv. minimální kostry v grafech. Zajímavostí je skutečnost, že existuje hned několik hladových algoritmů, které hledají minimální kostru grafu.
Mezi nejznámější patří algoritmus Borůvkův, Kruskalův a Jarníkův(Primův) (všimněte si, že dva z autorů jsou Češi). Základní přehled všech tří algoritmů podává česká Wikipedia.
Efektivita všech výše uvedených algoritmů, však závisí na implementaci elementárních kroků jako se setřídění hran, testování příslušnosti uzlů ke komponentám souvislosti resp. jejich spojování. Pro malé grafy lze využít triviální přístupy (označování uzlů a jejich pomalé prohledávání).
Knihovna NetworkX obsahuje implementaci Kruskalova algoritmu ve funkci minimum_spanning_tree resp. minimum_spanning_edges. Nalezenou minimální kostru je možné vizualizovat pomocí explicitního vykreslení hran
kostry tj. grafu, který vrací funkce minimum_spanning_tree.
Upozornění: Graf je náhodný a stejně tak jsou náhodné pozice uzlů získané algoritmem pro jeho vykreslení. Pokud není výsledek použitelný (hrany a uzly se příliš překrývají), zkuste kód vyhodnotit opakovaně.
from networkx.generators.random_graphs import fast_gnp_random_graph
from random import randint
import networkx as nx
import matplotlib.pyplot as plt
def plot_weighted_graph(g, mst):
pos = nx.spring_layout(g)
nx.draw_networkx_nodes(g, pos)
nx.draw_networkx_labels(g,pos)
nx.draw_networkx_edges(g, pos)
nx.draw_networkx_edges(mst, pos, width=2.0, edge_color="r")
labels = {(b,e) : g[b][e]["weight"] for b, e in g.edges()}
nx.draw_networkx_edge_labels(g, pos, labels)
plt.show()
g = fast_gnp_random_graph(10, 0.4)
for b,e in g.edges():
g[b][e]["weight"] = randint(0, 100)
plot_weighted_graph(g, nx.minimum_spanning_tree(g))
Úkol: Implementujte Borůvkův algoritmus, a vizualizujte jeho výsledek
Výše uvedená implementace není zdaleka optimální, a to ze dvou důvodů:
- opakovaně testuje hrany, které jsou již delší dobu součástí jediné komponenty
- algoritmus sjednocení (union) dvou komponent není efektivní
Úkol: Zkuste algoritmus zefektivnit. Rady: neoptimální hrany vedoucí k cyklům můžete odstranit, optimální datovou strukturou je union-find data structure (viz například https://www.geeksforgeeks.org/union-find)
Úkol: Implementujte Kruskalův algoritmus. Porovnejte jeho efektivitu s Borůvkovým pro různé grafy.
Hladové heuristiky¶
Hladové heuristické algoritmy, vycházejí z lokálních rozhodnutí, nemusí však najít skutečné globální optimum. I když ve většině případů najdou řešení, které je akceptovatelnou aproximací, existují úlohy a vstupy, pro něž určité heuristické algoritmy najdou jen špatná řešení (což jsou řešení, které jsou horší než náhodný přístup).
Heuristické algoritmy nacházejí uplatnění především v případech, kdy neexistuje řešení s polynomiální časovou složitostí, které je aplikovatelné jen na soubory s malým počtem dat. Typickým příkladem je problém obchodního cestujícího. Tento problém je NP-úplný, polynomiální algoritmus tudíž není znám a pravděpodobně ani neexistuje.
Základní heuristický algoritmus využívají i ti nejhloupější obchodní cestující: stačí vždy cestovat do nejbližšího nenavštíveného města (přičemž počáteční město je zvoleno náhodně, např. je to bydliště daného obchodního cestujícího). Tato heuristika vede v praktickém životě k akceptovatelným výsledkům (viz výsledek zdroje [22] na Wikipedii, který uvádí průměrné prodloužení o 25%).
Zajímavější heuristikou problémů obchodného cestujícího je využití optimalizace mravenčí kolonie, viz například implementaci na GitHubu.
Barvení grafu¶
Jedním z nejtypičtějších a nejpraktičtějších využití hladové heuristiky je barvení grafu. I když v případě https://mathworld.wolfram.com/VertexColoring.html existuje velké množství poznátků o optimálním barvení (optimální barvení uzlů využívá mimimální počet barev, což je tzv. chromatické číslo $\chi(G)$) neexistuje žádný efektivní algoritmus, který by pro libovolný graf našel optimální obarvení jinak než hrubou silou (přirozeně s exponenciální složitostí).
Nejčastěji se tak používá heuristika, která využivá hladový přístup. V nejjednodušším případě se stanoví nějaká posloupnost uzlů a v této posloupnosti se každý uzel obarví nejmenší možnou barvou (tj. předpokládá se že barvy jsou representovány přirozenými čísly).
Implementace této heuristiky je uvedena přímo v příslušném článku Wikipedie Greedy coloring.
Úkol: Ověřte různé strategie pořadí uzlů podporované funkcí zajišťující barvení hladovou heuristikou v knihovně NetworkX (viz https://networkx.github.io/documentation/stable/reference/algorithms/coloring.html). Vyzkoušejte pro náhodné grafy, či grafy speciální bipartitní, rovinné apod.
Úkol (k zamyšlení): Implementace funkce vracející nejmenší nepoužité číslo (tj. barvu nepoužitou sousedními uzly) je dosti neefektivní. V případě použití menších čísel (což lze u reálně použitelných grafů očekávat) lze využít bitovou aritmetiku. Zkuste najít příslusné řešení (rada: identifikace nejnižšího nenulového bytu a použití slovníku).
Další algoritmy využívající lokální rozhodování¶
Triviální řadící algoritmy¶
Jako triviální řadící algoritmy se označují řadící algoritmy, které mají časovou složitost $O(n^2)$ (tj. jejich časová složitost odpovídá porovnání každého prvku seznamu se všemi ostatními prvky).
Pro třídění větších polí nepředstavují konkurenci pro algoritmy založené na přístupu rozděl a panuj s časovou složitostí $O(n\log(n))$.
V praxi se však některé z nich mohou uplatnit ve specifických situacích:
malé kolekce (řádu malých desítek) = algoritmy insertion sort a v menší selection sort* jsou pro malé kolekce efektivnější než komplexní algoritmy (často se používají pro třídění listových úloh v přístupu rozděl a panuj)
snižení počtu výměn či zápisů = algoritmus selection sort a především cycle sort minimalizují počet záměn a tím i poče zápisů. To je výhodné,pokud je zápis výrazně pomalejší než čtení neb je počet zápisů omezen (flash paměti, EPROM paměti)
Úkol: Nastudujte výše uvedené řadící algoritmy a porovnejte jejich efektivitu pro malá pole.
Prohledávání v grafech¶
Klasickými algoritmy nad grafy využívajícími lokální volbu jsou algoritmy, které prohledávají všechny uzly dostupné z určitého uzlu (tj. v případě neorientovaného grafu celou komponentu).
Tuto třídu algoritmů lze využívat pro stromy i pro obecné grafy a lze je využívat pro různé účely od hledání nejkratších cest (= nejmenšího počtu procházených hran), přes zjišťování komponent po hledání různých charakteristik konektivity.
Dva základní algoritmy procházení grafů je prohledávání do šířky (breadth-firts) a do hloubky (depth-first).
Popis obou algoritmů můžete najít na https://www.educative.io/edpresso/how-to-implement-depth-first-search-in-python resp. https://www.educative.io/edpresso/how-to-implement-a-breadth-first-search-in-python.
Úkol: Implementujte algoritmus nad representací grafů použitou v knihovně NetworkX a porovnejte s implementací, která je již hotová v knihovně Networkx.
Úkol+: Rozšiřte algoritmus tak, aby rozlišil stromové, dopředné, zpětné a příčné hrany nad kostrou definovanou prohledáváním do hloubky (popis viz např. https://cw.fel.cvut.cz/old/_media/courses/a7b33tin/02_pruchod_grafem.pdf).
6. Algoritmy využívající náhodnost (pravděpodobnostní algoritmy)¶
Algoritmy využívající náhodnost se využívají v případech, kdy je dosažení určitého optimálního řešení neefektivní a postačuje jen řešení, které je pouze řešením přibližným (heuristické algoritmy), ale také v případech, kdy zahrnutí náhodnosti neovlivňuje výsledek, ale pouze efektivitu algoritmu. Algoritmus tak dosahuje „průměrné“ časové či paměťové složitosti bez ohledu na vstupní data. To je důležité především v případě, kdy je nejhorší časová složitost jiné třídy (průměr není jednoduchým aritmetickým průměrem) a kdy je pravděpodobnost výskytu kritických dat (= dat vedoucích k vyšší asymptotické složitosti) vysoká.
Příkladem elementární heuristiky je např. výpočet průměrů, který běžně vyžaduje součet všech prvků v kolekci (časová složitost je tedy $O(n)$). Pro zjištění odhadu průměru však stačí sečíst jen $k$ vybraných čísel (kde $k$ je konstanta $\ll n$), jehož časová složitost je konstantní. Je to stejný přístup, který se běžně používá ve statistickém šetření, které se provádí jen na malém vzorku populace (statistika nám při splnění určitých předpokladů dává i chyby takových odhadů).
Příkladem klasického pravděpodobnostního algoritmu, který snižuje pravděpodobnost extrémního času vykonávání je qsort s náhodnou volbou pivota, který dosahuje časové $O(log(n))$ i pro data seřazená v opačném gardu (tj. když např. chceme seřadit vzestupně data již seřazená sestupně). Pravděpodobnost takových data je výrazně vyšší než pravděpodobnost většiny ostatních permutací (která je při hodnotě $\frac{1}{n!}$ pro větší kolekce dat zanedbatelně malá). Důvodem je skutečnost, že data bývají často (částečně) setříděna/předtříděna v rámci předchozích zpracování (např. pokud vznikají spojením několika externích zdrojů).
Nejdříve se však podíváme na zdroj náhodných čísel tzv. generátory.
Hardwarové generátory náhodných čísel¶
Jádrem pravděpodobnostních algoritmů je získání náhodného proudu bitů s rovnoměrným rozdělením (tj. každý bit proudu má stejnou pravděpodobnost hodnoty 0 resp. 1) tj. například neexistuje žádná závislost mezi libovolnými dvěma bity (tj. například pravděpodobnost hodnoty v n-bitu není ovlivněna tím že $i*n$ bity pro $i\in{1..5}$ jsou jedničkové).
Získání skutečně náhodných posloupností je velmi obtížné. I když je podle našeho současného fyzikálního modelu produkují některé kvantové procesy, je jejich získávání technicky dosti obtížné (včetně nutnosti transformací, pokud nemají rovnoměrné rozložení). Prakticky se proto používá šum, který vzniká překrýváním většího počtu kvantových procesů (termální šum, fotoelektrický efekt a podobně).
Reálný generátor využívající radioaktivníh rozpad cesia $^{137}\mathrm{Cs}$ můžete najít na adrese https://www.fourmilab.ch/hotbits/.
V případě pravděpodobnostních algoritmů však hardwarové zdroje nehrají téměř žádnou úlohu, neboť jsou:
- příliš pomalé (tj. čas produkce náhodných čísel nelze zanedbat, výše uvedený generátor například
hotbitsprodukuje jen cca 100 bytů/s) - nejsou reprodukovatelné (a tudíž obtížněji laditelné).
Hardwarové generátory náhodného bitového proudu je však nutno použít tehdy, pokud hrozí nebezpečí, že zpracovávaná data mohou být vědomě ovlivněna, tak aby s predikovanými reprodukovatelnými "náhodnými" hodnotami pravděpodobnostní hodnoty vracely chybné výsledky resp. spotřebovával enormně mnoho prostředků.
Útočník může například přeuspořádat vstupní pole tak, že operace shuffle (proházení prvků) nevrátí (zdánlivě) náhodný mix prvků, ale seřazené pole. Pokud na toto pole použijeme quicksort je to nejen nadbytečné, ale může vést k tomu že časová složitost algoritmu může dosáhnout $O(n^2)$, což může tuto operaci výrazně zpomalit. I když se útoky tohoto typu mohou jevit v této oblasti jako nerealistické, nelze je vyloučit v případě, kdy suboptimální výkonnost může snížit dostupnost nějaké služby. V tomto případě je nutno na reprodukovatelnost rezignovat (klíčová je jen ve fázi ladění) a použít hybridní přístup (hardwarově generovaný proud slouží pro občasnou inicializaci softwarových generátorů).
V operačním systému Unix lze pro tyto účely využít objekt třídy random.SystemRandom, který spojuje softwarový generátor se zdroji tzv. entropie (= zdroje šumu + hardwarové generátory, které jsou dnes k dispozici i v rámci běžných CPU).
Upozornění: I když tento generátor není reálně reprodukovatelný, jeho použití pro kryptografické účely není doporučováno (Python poskytuje pro kryptografické účely modul
secrets)
from random import SystemRandom
generator = SystemRandom()
print(bin(generator.getrandbits(1024))) # využívá systémový generátor náhodných bitů pro získání produ 1024 bitů (= 128 bytů)
Transformace na posloupnost čísel¶
Než přejdeme k softwarovým generátorům je nutno prodiskutovat mapování proudu bitů na proud náhodných čísel a to jak celých (typu int) tak čísel pohyblivou řádovou čárkou (typ float). Nejjednodušší situace je v případě celých čísel v rozsahu 0 (včetně) až $2^n$ (vyjma). Zde stačí vzít $n$ bitů a interpretovat je jako číslo. V případě, že výsledkem mají být celá čísla v jiném rozsahu, je nutný přepočet tj. například aplikace operace zbytku po dělení (pokud má být výsledkem celé číslo) resp. dělení. Ti však má ve většině případů bohužel tyto negativní důsledky:
- výsledné rozdělení čísel nemusí být rovnoměrné
- při transformaci se zahazují určité bity (tento nedostatek nemá vliv na výsledky pravděpodobnostních algoritmů, snižuje však jejich efektivitu)
- používají neefektivní operace (například transformaci mezi celými čísly a čísly v řádové čárce)
Elementární přístupy jako je vedou vždy k alespoň jednomu z výše uvedených negativních důsledků a často i ke dvěma.
Úkol: Navrhněte triviální a implemntujte algoritmus pro generování celých čísel v rozsahu 0 až $n$, na základě možnosti generovat libovolnou posloupnost náhodných bitů (viz výše metoda
getrandbits). Diskutujte řešení.
Nejjednoduší přístup je výpočet podle vztahu i % n, kde i je náhodné číslo v rozsahu 0 az $2^k$, přičemž $2^k > n$.
Je-li $2^k = n$: triviální případ, řešení je přesné a rychlé i úsporné (není potřeba ani provádět zbytek po dělení)
Je-li $2^k$ řádově shodné s $n$: například $n=6$ a $k$ je 3, 4, apod., pak je výpočet efektivní, ztrácí se jen pár bitů ale výsledné rozložení pravděpodobnosti není rozhodně rovnoměrné. Při volbě $k=3$ budou mít čísla 0 a 1 dvakrát vyšší pravděpodobnost než čísla 2 až 5 (tj. simulovaná kostka bude viditelně nepřející).
Je-li $2^k \ll n$:
odchylky od rovnoměrného rozložení pravděpodobnosti bude malé (a ve většině případů zanedbatelné), mnoho bitů však přijde nazmar. Navíc výpočet 2**k % n může vyžadovat více procesorového času pro větší $k$.
Alternativně lze využít vypočet prostřednictvím typu float v podobě int(i/(2**k)*n, která trpí stejnými vadami (např. pro nízké hodnoty 2^k je rozdělení podobně nepravidelné, i když méně zřejmé). Využití dělení v plovoucí řádové čárce a přetypování na int může být značně neefektivní (především pokud cílová platforma postrádá efektivní hardwarové FPU.
Úkol: vypočítejte rozdělení pravděpodobnosti hodnot 0-5 vypočítaných pomocí vztahu, kde $k=3,8,16,32$ (tj. nejmenší možné a o délce slov běžných procesorů). Zamyslete se nad výsledkem (jak nepravidelnosti tolerovatelné, resp, v jakých kontextech)
from collections import Counter
from matplotlib.pyplot import bar, show
for k in [3,8,16]:
nk = 2**k
c = Counter(int(i/nk*6) for i in range(nk))
h = [(n, 100*p/nk) for n,p in c.most_common()]
x,y = zip(*h)
print(k, h)
bar(x,y)
show()
Softwarové generátory pseudonáhodných čísel¶
Sigtwarové generátory produkují náhodná čísla tím, že interně udržují nějaký stav (pole bytových hodnot), na který je možno aplikovat funkci $t_s$, která vrátí nový stav.
$S_{i+1} = t_s(S_i)$, kde $i \geq 0$
Počáteční hodnota $S_0$ se se získává transformací ($t_i$) z hodnoty nazývané inicializační (random seed), která může být buď fixní (umožňujíc reprodukovatelnost), měnit se s časem (typicky je to např. některá z representací reálného času) nebo to může být náhodná hodnota získaná hardwarovým generátorem.
Aplikací další funkce $t_r$ na jednotlivé stavy pak získáváme postupně bloky pseudonáhodných bitů (které jsou typicky chápány jako čísla)
$r_i = t_r(S_i)$
Ze zjednodušeného předpisu je zřejmé, že výsledkem aplikace softwarového generátoru nemůže být skutečně náhodná posloupnost bitů resp. čísel. Neboť hodnoty $r_i$ jsou evidentně závislé na předchozích stavech resp. na inicializační hodnotě (lze je tudíš snadno znovu reprodukovat). Navíc se hodnoty musí cyklicky opakovat (stav je representovaný konečým polem bitů a tak existuje jen konečný počet stavů).
Na druhou stranu nelze v mnoha případech (a v konečném a rozumném čase rozlišit) mezi posloupností genrovanou hardwarovým generátorem a posloupností generovanou generátorem softwarovým, neboť:
- perioda po níž se čísla cyklicky opakují je mimořádně dlouhá
- hodnoty jsou rovnoměrně rozloženy (tj. nelze statisticky dokázat, že existuje nějaký bias)
- libovolné dostatečně velké nepřekrývající se výběry nejsou korelované (to platí i pro hodnoty, které následují bezprostředně či blízko za sebou)
- vzdálenosti mezi opakujícími se vzory mají stejnou distribuci jako v případě skutečně náhodných posloupností
- hodnoty jsou rovnoměrně distribuovány i ve více dimenzích (tj. body s s pseudonáhodně gnerovanými souřadnicemi vyplňují rovnoměrně každý či alespoň daný $n$-dimezionální prostor).
Lineární kongruentní generátory¶
Nejjednoduššími softwarovými generátory jsou lineární kongruenční generátory, u nichž je následující číslo generováno předpisem:
$x_{i+1}=(ax_{i}+c)\ {\bmod {\ }}m$, kde $a$, $c$ a $m$ jsou vhodně zvolená přirozená čísla.
I když jsou známy základní požadavky na volbu čísel, které mimo jiné zajistí, že cyklus bude mít maximální délku $m-1$, je nalezení skutečně optimálních čísel spíše empirickou úlohou. Navíc je výhodné, pokud $m$ odpovídá šířce slova (tj. je rovna např. 2^32 u 32-bitových procesorů), neboť to usnadňuje implementaci zbytku po dělení (operace modulo je eliminována neboť všechny výpočty nativně probíhají v modulární aritmetice $2^n$, tj. vyšší bity výsledku jsou automaticky ignorovány) a/nebo $c$ c je rovno nule (eliminace sčítání) resp. 1 (inkrementace je na některých architektur efektivnější než obecný součet). Lineárně kongruentní generátory jsou proto obecně velmi rychlé (tím spíše pokud architektura podporuje rychlou operaci multiply-accumulation).
Poznámka: Některé volby parametrů vedou k tomu, že některé bity mají nižší periodicitu a tak je nutno ze stavu extrahovat jen některé bity (typické především tehdy je-li parametr $c$ nulový).
Úkol: implementujte tři typické lineární kongruenční operátory a diskutujte jejich použitelnost.
- klasická implementace operace random v jazyce C: a = 1103515245, c = 12345, m = 2^31 (typicky se používá bity vyšší bity (16-30).
- minimalistická implementace (CRAY RANF): a = 44 485 709 377 909, c = 0, m = 2^48 - 1
- katastrofická implementace (RANDU): a = 65539, c = 0, m = 2^31
Otázka pro zvídavé: Aditivní člen v prvním příkladě má hodnotu 12345. Je to z důvodu, že implementátory jiné „náhodné“ číslo nenapadlo?
from itertools import islice
def lcg(modulus, a, c, seed):
"""Linear congruential generator.
source: https://en.wikipedia.org/wiki/Linear_congruential_generator"""
while True:
seed = (a * seed + c) % modulus
yield seed
for number in islice(lcg(2**31, 1103515245, 12345, 0), 12):
print(f"{number:032b}")
Linearní kongruentní generátory jsou vhodné, jen tehdy pokud:
- implementace musí využívat jen velmi omezené prostředky (např. na známých Arduinech s pár KiB paměti a pomalými aritmetickými operacemi)
- postačuje menší počet generovaných čísel (což často neplatí u Monte Carlo algoritmů!), u 32 bitové implementace je vhodný například pro řády tisíců náhodných hodnot
- jsou využívány jen vyšší bity stavu (což ale snižuje jejich efektivitu). Zvlášť problematický je nejnižší bit (viz příklad výše)
Mersenne twister a PCG gnerátory¶
Nejpopulárnějšm algoritmem dneška je tzv. Mersenne twister, který využívá mnohem větší representaci stavů a řeší většinu problémů LCG generátorů. Zvláště efektní je extrémně dlouhá perioda, která u běžně používané konfigurace dosahuje $2^{19937}-1 \approx 10^6000$, což by mělo stačit po celou zbývající dobu existence lidské civilizace.
Na druhou stranu má více parametrů, které je nutno pečlivě volit, vyžaduje více paměti a je pomalejší (i když nikoliv řádově).
Mersenne twister je standardním softwarovým generátorem, používaným v pythonském modulu random (platné pro Python 3.8 a starší verze Pythonu 3) a to hned k několika účelům:
- generování posloupnosti náhodných bitů
random.getrandbits - generování náhodných celých čísel z určitého rozsahu (
random.randintresp.random.randrange) - náhodný výběr ze seznamu (
random.choicesneborandom.sampleresp. specializovanérandom.choice) - náhodná permutace prvků seznamu (
shuffle) - generování čísel v pohyblivé řádové s rovnoměrným rozdělením čárce z intervalu $[0, 1)$ (
random.random) resp. z intervalu $[a, b]$ (random.uniform) - generování čísel v jinách spojitých rozděleních (včetně normálního a ecponenciálního, viz níže)
Jediné, co funkce z modulu random přímo nenabízí je generování kolekcí hodnot (to lze samozřejmě zajistit opakovaným voláním, ale to je neefektivní). V tomto případě je vhodnější využít modul numpy.random.
Při použití tohoto modulu je možno (od verze 1.18) využít dvou přístupů:
- starší (klasický) přístup (jednoduché volání funkcí, jejichž jediným parametrem je tvar výsledného pole), modul nabizí velké množství funkcí pro různé typy výstupů, různá rozdělení (a to v často v různých variantách, aby bylo dosaženo kompatibility s jiným matematickým softwarem). Základním generátorem je klasický Mersenne twister.
- nový přístup (oddělení bitového generátoru a generátoru čísel), přehlednější a ortogonálnější rozhraní, explicitní podpora pro vytváření nezávislých generátorů v paralelních aplikacích (tzv. spawning). Navíc nabízí jako standardní základní (bitový) generátor relativně nový PCG64 generátor, který by měl mít lepší statistické vlastnosti při využití méně prostředků (včetně času). Pro ty, kteří novinkám nevěří, je k dispozici standardní Mersenne Twister.
Popis PCG generátorů naleznete na stránce https://www.pcg-random.org/index.html. Obsahují i zajímavý úvod ke statistickému testování generátorů a také diskutují silné i slabé stránky většiny v současnosti využívaných generátorů.
Úkol:
- vyplňte seznam 100 náhodnými hodnotami typu
intz intervalu 1 .. 100 (s využitím jediného příkazu a klasické knihovnyrandom)- vyplňte matici 100×100 náhodnými hodnotami typu
floatz intervalu $[-100, 100]$ (s využitím starého i nového rozhraní v NumPy).
Generování náhodných hodnot z jiných rozdělení¶
Častým požadavkem je generování pseudonáhodných čísel z jiného než rovnoměrného rozdělení (typicky jsou to spojitá rozdělení a generována jsou tak čísla typu float). Pro tyto účely existují dva přístupy, které se typicky kombinují. Oba využívají jako zdroj generátor s náhodným rozdělením vracející čísla float z určitého intervalu (standardně $[0, 1)$).
- na hodnoty je applikována nějaká (nelineární) funkce
- některé hodnoty jsou zahozeny, tak aby měl výsledný generátor příslušné rozdělení
V případě rozumných generátorů je transformační funkce jednoduchá (a tudíž rychlá) a počet zahozených hodnot malý (v opačném případě by algoritmus produkoval jen málo čísel za časovou jednotku a zpomaloval by tak příslušný prevděpodobnostní algoritmus).
V minulosti se používalo velké množství různých algoritmů. Dnes se však využívá tzv. zikkurat algoritmus, pojmenovaný podle mezopotámských městských chrámů.
Přehledný a přitom dostatečně detailní popis algoritmu (pro normální rozdělení) najdete na https://heliosphan.org/zigguratalgorithm/zigguratalgorithm.html (naleznete tam i obrázek, který vysvětluje název algoritmu). Algoritmus je elegantní a jednoduchý vyjma speciálního ošetření tzv. ocasu (část křivky hustoty rozdělení, kde je hodnota hustoty pravděpodobnosti nízká, ale nikdy nedosáhne nuly).
Úkol: Implementujte zikkurat algoritmus pro symetrické trojúhelníkové rozdělení (tj. modus leží uprostřed mezi krajními body, vně jichž je hustota pravděpodobnosti nulová).
Úkol: Vygenerujte 1000 bodů s exponenciálním rozdělení se střední hodnotou 0.5 a zobrazte jejich histogram (s krokem 0.1) a proložte ho teoretickou křikvkou hustoty
Poznámka a úkol: Výše uvedený přístup lze použít i obráceně. Mějme generátor binárních hodnot, který má tzv. bias tj. častěji generuje jedno ze dvou binárních čísel častějí (tj. s pravděpodobností $p \neq 0.5$). Implementujte řešení ze stránky https://www.geeksforgeeks.org/print-0-and-1-with-50-probability/ a diskutujte ho.
Randomizované verze standardních algoritmů¶
Jedním z využití generátoru pseudonáhodných čísel je tzv. randomizované verze algoritmů, tj. verze v nichž dochází buď k randomizaci vstupu resp. k randomizaci nějakého kroku algoritmu.
Prvním příkladem je randomizace vstupu při vytváření binárního stromu. Tato strategie je vhodná v případě, že ukládáno je velké množství dat najednou a přidávání jednotlivých dat je jen výjimečné.
Algoritmus se v tomto případě nemění, vstupní data však musí být dočasně uložena v poli či seznamu a náhodně promíchána (shiffling).
Nejčastěji používaným algoritmem používaným pro náhodné promíchání je algoritmus Fisher-Yattesův, který nabízí promíchání s časovou složitostí $O(n)$ a probíhá přímo ve zdrojovém poli {in-place).
Úkol: Prostudujte moderní verzi Fisher-Yettesova algoritmu na https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle (velmi kvalitní Wikipedia článek, diskutující i jiné přístupy a také úskalí, které může být dáno např. použitím nedokonalého generátoru náhodných čísel). Algoritmus efektivně implementujte.
Fisher-Yattesův algoritmus je přirozeně implementován i v ve standardní knihovně Pythonu (random.shuffle) i v Numpy pro pole (numpy.random.shuffle).
Úkol k zamyšlení: V případě jednorázové inicializace vyhledávacího stromu (tj. jedno vkládání, mnoho vyhledávání) lze alternativně využít setřídění a následné využití binárního vyhledávání. Diskutujte prostorovou a časovou složitost.
Příkladem interně ranodmizovaného algoritmu je quicksort s náhodným výberem pivota. Stránka https://www.geeksforgeeks.org/quicksort-using-random-pivoting/ obsahuje dvě různé implementace využívající
různá schémata dělení (porovnejte!).