10: Paralelní výpočty a profilování
Tento notebook je výukovým materiálem v předmětu BI-JUL.21 vyučovaném v zimním semestru akademického roku 2021/2022 Tomášem Kalvodou. Tvorba těchto materiálů byla podpořena NVS FIT.
Hlavní stránkou předmětu, kde jsou i další notebooky a zajímavé informace, je jeho Course Pages stránka.
versioninfo()
1. Úvod
V Julia prostředí se můžeme setkat s čtyřmi druhy paralelismu:
- asynchroní "tasky" (coroutines),
- více vláknové výpočty (multi-threading),
- distribuované výpočty (distributed computing),
- výpočty na GPU.
V tomto notebooku se zaměříme na druhý bod a třetí bod. Prezentace bude více experimentální. Zájemci mohou prozkoumat následující balíčky a zdroje
2. Paralelní výpočty
2.1: Více vláknové výpočty
Julia nám umožňuje mezi vlákna rozdělovat tasky, vlákna sdílí paměťové prostředky.
Z příkazové řádky můžeme Julia spustit s více vlákny pomocí přepínače -t
/--threads
:
$ julia -t4
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.6.3 (2021-09-23)
_/ |\__'_|_|_|\__'_| |
|__/ |
julia> Threads.nthreads()
4
Pokud chceme pracovat s více vlákny i v Jupyter notebooku/labu tak musíme nastavit proměnnou prostředí JULIA_NUM_THREADS
na požadovaný počet vláken.
Pokud toto provedete (například pomocí export JULIA_NUM_THREADS=4
na Linuxu), pak by vám výstup následující buňky měl dát 4.
Threads.nthreads()
Pokud jste spustili Jupyter notebook/lab standardním způsobem, pravděpodobně je výstup předešlé buňky 1.
Threads.@threads
makro
Toto makro z modulu Threads
slouží k jednoduchému paralelizování for
cyklu.
Pokud ho aplikujeme na for
cyklus, tak Julia rozdělí práci rovnoměrně (co do počtu prvků v iterátoru) mezi jednotlivá vlákna a počká až všechny vykonají svou práci.
Aktuálně (tj. zima 2021) je k dispozici pouze toto "statické" rozdělení práce, v budoucnou přibudou i další možné strategie, což se už stalo (viz dokumentaci zmíněného makra, k dispozici máme přepínače :dynamic
a :static
).
Toto je samozřejmě vhodné pouze pokud každý z úkonů v cyklu trvá přibližně stejně dlouho.
Jako jednoduchý názorný příklad uvažme pole nul a do každého zapišme id vlákna. Nejprve znovu zkontrolujme, že nám běží jádro s více vlákny:
Threads.nthreads()
Julia vlákna interně čísluje od 1. Demonstrační příklad je tedy následující:
zrs = zeros(Int64, 10)
for j in 1:10
zrs[j] = Threads.threadid()
end
zrs
# pole o deseti prvcích
L = 10
zrs = zeros(Int64, L)
Threads.@threads :dynamic for j in 1:L
# do každého prvku zapíšeme id vlákna
zrs[j] = Threads.threadid()
end
zrs
# pole o deseti prvcích
L = 10
zrs = zeros(Int64, L)
Threads.@threads :static for j in 1:L
# do každého prvku zapíšeme id vlákna
zrs[j] = Threads.threadid()
end
zrs
Co když každý krok trvá jinou dobu?
# pole o deseti prvcích
L = 20
zrs = zeros(Int64, L)
Threads.@threads :dynamic for j in 1:L
sleep(rand()/2)
# do každého prvku zapíšeme id vlákna
zrs[j] = Threads.threadid()
end
zrs
# pole o deseti prvcích
L = 20
zrs = zeros(Int64, L)
Threads.@threads :static for j in 1:L
sleep(rand()/2)
# do každého prvku zapíšeme id vlákna
zrs[j] = Threads.threadid()
end
zrs
Threads.@spawn
makro, fetch
a wait
@spawn
makro z modulu Threads
dokáže zadat práci nějakému vláknu bez čekání na výsledek.
Jeho návratová hodnota je typu Task
, pomocí které můžeme později výsledek získat pomocí metody fetch
nebo počkat až výpočet doběhne pomocí metody wait
.
using BenchmarkTools
Threads.nthreads()
Jako demonstrační příklad si představme následující ospalou funkci:
"""
Chvilku spinkej a pak vrať 42!
"""
function sleepy_function(t=2)
sleep(t)
return 42
end
sleepy_function()
Představme si, že chceme zavolat tuto funkci tolikrát, kolik máme k dispozici vláken a poté výsledky sečíst. Nejprve sériově:
function f()
values = zeros(Int64, Threads.nthreads())
for j = 1:Threads.nthreads()
values[j] = sleepy_function()
end
return sum(values)
end
Jaký je výsledek?
f()
Kolik času to zabralo?
@benchmark f()
Osm sekund je přesně to co jsme čekali. Pojďme tuto funkci nyní skutečně paralelizovat a to rovnou dvěma způsoby:
"""
@spawn-fetch verze.
"""
function g()
values = Array{Task}(undef, Threads.nthreads())
for j = 1:Threads.nthreads()
values[j] = Threads.@spawn sleepy_function()
end
return sum(fetch.(values))
end
"""
@threads verze (:dynamic).
"""
function h()
values = zeros(Int64, Threads.nthreads())
Threads.@threads :dynamic for j = 1:Threads.nthreads()
values[j] = sleepy_function()
end
return sum(values)
end
"""
@threads verze (:static).
"""
function ch()
values = zeros(Int64, Threads.nthreads())
Threads.@threads :static for j = 1:Threads.nthreads()
values[j] = sleepy_function()
end
return sum(values)
end
g()
h()
ch()
@benchmark g()
@benchmark h()
@benchmark ch()
Opět, dvě sekundy na čtyřech vláknech jsou přesně to co jsme čekali.
Cvičení, resp. poznámka k 1. domácímu úkolu (z roku 2021)
Porovnejme rychlost několika způsobů jak k jednomu řádku matice přičíst násobek jiného řádku matice (úloha z prvního domácího úkolu). Rozebrme rozdíly mezi jednotlivými přístupy.
function add1!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
m[target, :] += α * m[source, :]
end
function add2!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
m[target, :] += α .* m[source, :]
end
function add3!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
for j in axes(m, 2)
m[target, j] += α * m[source, j]
end
end
function add4!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
@inbounds for j in axes(m, 2)
m[target, j] += α * m[source, j]
end
end
function add5!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
@simd for j in axes(m, 2)
@inbounds m[target, j] += α * m[source, j]
end
end
function add6!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
@simd for j in axes(m, 2)
m[target, j] += α * m[source, j]
end
end
function add7!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
Threads.@threads :dynamic for j in axes(m, 2)
@inbounds m[target, j] += α * m[source, j]
end
end
function add8!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
Threads.@threads :static for j in axes(m, 2)
@inbounds m[target, j] += α * m[source, j]
end
end
using BenchmarkTools
mat = 20 * (rand(4, 10^8) .- 1);
@benchmark add1!(a, -2.376, 1, 2) setup=(a=copy(mat))
@benchmark add2!(a, -2.376, 1, 2) setup=(a=copy(mat))
@benchmark add3!(a, -2.376, 1, 2) setup=(a=copy(mat))
@benchmark add4!(a, -2.376, 1, 2) setup=(a=copy(mat))
@benchmark add5!(a, -2.376, 1, 2) setup=(a=copy(mat))
@benchmark add6!(a, -2.376, 1, 2) setup=(a=copy(mat))
@benchmark add7!(a, -2.376, 1, 2) setup=(a=copy(mat))
@benchmark add8!(a, -2.376, 1, 2) setup=(a=copy(mat))
Jaký učiníme závěr z těchto experimentů?
Cvičení: Paralelizace výpočtu Julia fraktálu
Pojďme se vrátit k dřívějšímu příkladu na výpočet Julia fraktálu (notebook bi-jul-07.ipynb
).
Aby se nám kód dobře paralelizoval, tak ho lehce rozdělíme:
using PyPlot, BenchmarkTools
"""
Výpočet počtu iterací `f` v bodě `z` nutných k opuštění koule o poloměru `R` a středu
v počátku komplexní roviny.
"""
function julia_iterations(z::Complex{Float64}, f::Function, R::Float64; imax::Int64=1000)
for n = 1:imax
z = f(z)
if abs(z) > R
return n
end
end
return imax
end
"""
Sériový výpočet Julia fraktálu.
"""
function julia_matrix!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
f(u) = u^2 + c
for j = axes(reals, 1), k = axes(imags, 1)
data[j, k] = julia_iterations(reals[j] + 1im * imags[k], f, R, imax = imax)
end
end
Následující výpočet zatím není paralelizovaný, ověřme jeho základní funkčnost.
res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix!(-0.8 + 0.256im, 10.0, res, ims, data)
fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])
A podívejme se na dobu běhu.
@benchmark julia_matrix!(-0.8 + 0.256im, 10.0, $res, $ims, $data)
Kolik máme k dispozici vláken? Případně nastavte proměnnou prostředí JULIA_NUM_THREADS
před spuštěním Jupyter notebooku/labu.
Threads.nthreads()
Paralelizace se v tomto případě zdá poměrně přímočarou.
Ovšem každý "job" (iterace ) nepoběží úplně stejně dlouho.
Pozor, Threads.@threads
zatím nepodporuje "složené cykly", musíme kód přepsat pomocí dvou vnořených for
cyklů (tj. paralelizujeme plnění matice po řádcích):
function julia_matrix_p1!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
f(u) = u^2 + c
Threads.@threads for j = axes(reals, 1)
for k = axes(imags, 1)
data[j, k] = julia_iterations(reals[j] + 1im * imags[k], f, R, imax = imax)
end
end
end
res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p1!(-0.8 + 0.256im, 10.0, res, ims, data)
fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])
@benchmark julia_matrix_p1!(-0.8 + 0.256im, 10.0, $res, $ims, $data)
Získali jsme přibližně trojnásobné zrychlení.
function julia_matrix_p2!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
f(u) = u^2 + c
Threads.@threads :static for j = axes(reals, 1)
for k = axes(imags, 1)
data[j, k] = julia_iterations(reals[j] + 1im * imags[k], f, R, imax = imax)
end
end
end
res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p2!(-0.8 + 0.256im, 10.0, res, ims, data)
fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])
@benchmark julia_matrix_p2!(-0.8 + 0.256im, 10.0, $res, $ims, $data)
Zajímavé, zkusme změnit parametry.
@benchmark julia_matrix_p1!(-0.8 + 0.256im, 5.0, $res, $ims, $data, imax=5_000)
@benchmark julia_matrix_p2!(-0.8 + 0.256im, 5.0, $res, $ims, $data, imax=5_000)
Zkusme aleternativní přístup pomocí @sync
a Threads.@spawn
. V tomto případě opět ručně iterujeme přes řádky matice.
function julia_matrix_p3!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
f(u) = u^2 + c
nthreads = Threads.nthreads()
@sync for t in 1:nthreads
Threads.@spawn for j = t:nthreads:length(reals), k = axes(imags, 1)
data[j, k] = julia_iterations(reals[j] + 1im * imags[k], f, R, imax = imax)
end
end
end
res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p3!(-0.8 + 0.256im, 10.0, res, ims, data)
fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])
@benchmark julia_matrix_p3!(-0.8 + 0.256im, 5.0, $res, $ims, $data, imax=5_000)
To je prakticky stejný výsledek, o fous rychlejší!
res = LinRange(-0.4, 0.4, 2000)
ims = LinRange(-0.6, -0.1, 1600)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p3!(-0.8 + 0.156im, 300.0, res, ims, data, imax=1000)
fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])
res = LinRange(-0.4, 0.4, 2000)
ims = LinRange(-0.6, -0.1, 1600)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p1!(-0.8 + 0.156im, 300.0, res, ims, data, imax=1000)
fig, ax = plt.subplots()
ax.imshow(transpose(data), extent=[res[1], res[end], ims[1], ims[end]])
2.2 Distribuované výpočty
Pokud chceme využít k výpočtům více strojů, například celé clustery, můžeme k tomu použít modul Distributed
ze standardní knihovny.
Julia procesy pak běží odděleně (potenciálně i na jiných strojích) a kontrola běhu programu se stává obtížnější/zábavnější.
Procesy jsou opět číslovány od jedné. Ten s id rovným jedné odpovídá procesu v němž běží REPL/Jupyter (interaktivní sezení).
Pokud chceme spustit Julia s více procesy (workery) na jednom stroji, použijeme k tomu přepínač -p
/--procs
.
Například požadujeme-li čtyři procesy na jednom stroji:
$ julia -p 4
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.6.3 (2021-09-23)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> nworkers()
4
Pokud chceme spustit více procesů na více strojích, musíme tyto uvést v souboru a ten předat Julia při spuštění pomocí argumentu --machine-file FILE
. Pro komunikaci se používá ssh protokol (bez hesla, klíč), řádek v daném souboru má formát [count*][user@]host[:port] [bind_addr[:port]].user
(výchozí hodnota user
je aktuální uživatel).
Alternativně lze procesy přidávat/odebírat programově pomocí metod addprocs
/rmprocs
z modulu Distributed
.
Viz jejich dokumentaci.
Pokud Julia spustíme tímto způsobem, tak automaticky dojde k importu modulu Distributed
.
Logika kontroly výpočtu je podobná jako u více vláknových výpočtů. Máme k dispozici následující makra/metody:
nworkers()
: počet procesů (workerů),workers()
: pole s id workerů,@spawnat ID EXPRESSION
: vyhodnotí výrazEXPRESSION
na procesu s idID
(může být rovno:any
, pokud nám nezáleží na kterém), vrací objekt typuFuture
.fetch(future)
,wait(future)
: získá výsledek vzdáleného výpočtufuture
, resp. počká na doběhnutí výpočtu.@everywhere EXPRESSION
: vyhodnotí výrazEXPRESSION
na všech procesech (například funkce, s kterými pracujeme, musí být dostupné všem procesům).pmap(f, array)
: paralelní map.@distributed [OP] for
: paralelní for cyklus s případnou redukcí operátoremOP
.
V tento okamžik doporučuji restartovat jupyter notebooku a spustit ho s jedním vláknem.
Threads.nthreads()
Demonstrační ukázka: náhodná procházka
Aktuálně nám v Jupyter notebooku pravděpodobně běží Julia s jedním procesem/workerem:
using Distributed
nworkers()
Threads.nthreads()
Zvyšme počet procesů na 4:
addprocs(4)
nworkers()
K demonstraci paralelizace zvolme modelování tzv. diskrétní náhodné procházky. Představme si opilého chodce, který vyráží na procházku z bodu na reálné ose a v každém kroku si s pravděpodobností 0.5 vybere jednotkový krok vpravo, nebo vlevo.
Do jaké vzdálenosti od počátku se dostane po krocích? Jak bude vypadat hustota pravděpodobnosti, že se po těchto krocích nachází v jistém bodě ?
K experimentálnímu zodpovězení těchto otázek budeme potřebovat tyto procházky simulovat mnohokrát. To je ideální příležitost pro paralelizaci. Sestrojme nejprve sériové řešení:
function random_walk(steps)
position = 0
for _ in 1:steps
if rand() >= 0.5
position += 1
else
position -= 1
end
end
return position
end
random_walk(0)
random_walk(10)
Vygenerujme sériově 10 takovýchto náhodných procházek o 1000 krocích.
function get_walks(n, steps)
return [ random_walk(steps) for _ in 1:n ]
end
get_walks(10, 1_000)
Kolik času nám zabere procházek o stejném počtu kroků?
@time get_walks(10_000, 1_000_000);
Nejprve musíme naší funkci vykonávající náhodnou procházku zadefinovat i v ostatních procesech. Teď by v našem procesu výpočet prošel, ale v jiném ne:
t = @spawnat 1 random_walk(100) # v tento okamžik OK
fetch(t)
t = @spawnat 2 random_walk(100) # v tento okamžik selže
fetch(t)
Stačí použít makro @everywhere
(typicky můžeme udělat třeba @everywhere include("common.jl")
, kde v daném souboru máme kód pořebný k provádění výpočtů).
@everywhere function random_walk(steps)
position = 0
for _ in 1:steps
if rand() >= 0.5
position += 1
else
position -= 1
end
end
return position
end
# @everywhere include("source.jl")
Pro jistotu zkontrolujme, že vše funguje tak jak očekáváme.
t = @spawnat 3 random_walk(100)
fetch(t)
K samotné paralelizaci použijeme metodu pmap
nebo paralelní for
cyklus (makro @distributed
).
pmap
rozděluje práci prvek po prvku a hodí se proto na úlohy, které jsou větší, nebo nerovnoměrné.
Naopak @distributed for
rozdělí práci mezi workery najednou.
test = zeros(10)
pmap(x -> (sleep(0.2); myid()), test)
workers()
myid()
function get_walks_p1(n, steps)
output = fill(steps, n)
output = pmap(random_walk, output)
return output
end
get_walks_p1(10, 1_000)
Kolik času nám zabere procházek o stejném počtu kroků?
@time get_walks_p1(10_000, 1_000_000);
Při použití paralelního for cyklu potřebujeme mít možnost zaznamenávat výsledky do sdílené paměti.
Přesně k tomu slouží SharedArrays
.
using SharedArrays
function get_walks_p2(n, steps)
output = SharedArray{Int64}(n)
@sync @distributed for j = 1:n
output[j] = random_walk(steps)
end
return output
end
get_walks_p2(10, 1_000)
Kolik času nám zabere procházek o stejném počtu kroků?
@time get_walks_p2(10_000, 1_000_000);
Proč byl bylo řešení s pmap
nejpomalejší? Protože úloha, kterou řešíme byla příliš jednoduchá a výhoda paralelizace byla utopena v komunikaci mezi procesy (pmap
distribuje každou z úloh zvlášť).
Ukažme si tento efekt na ještě jednodušším příkladě, uvažme jednoduchý výpočet kvadrátu a faktorizaci.
@everywhere begin
using Primes
f(n) = n^2 # "levná" operace
g(n) = factor(n) # občas "drahá" operace
end
A tři různé způsoby výpočtu těchto funkcí v jistém rozsahu.
function serial(func, rng)
map(func, rng)
end
function para2(func, rng)
pmap(func, rng)
end
function para3(func, rng)
@sync @distributed for j in rng
func(j)
end
end
integer_range = BigInt(10^18):BigInt(10^18 + 1_000)
Výpočty kvadrátů:
@time serial(f, integer_range);
@time para2(f, integer_range);
@time para3(f, integer_range);
Výpočty faktorizací.
@time serial(g, integer_range);
@time para2(g, integer_range);
@time para3(g, integer_range);
Vraťme se zpět k naší náhodné procházce, napočtětem 10000 procházek o milionu krocích (hodně bloudící opilec):
ws = get_walks_p2(10_000, 1_000_000)
using PyPlot
hist(ws, bins=1000);
Chudák opilec se tedy bude nejpravděpodobněji nacházet v bodě 0. Spočtěme si pár statistických ukazatelů.
using Statistics
mean(ws)
median(ws)
std(ws)
minimum(ws), maximum(ws)
Cvičení: Collatz conjecture
Collatzova hypotéza je pravděpodobně jedním z nejznámějších nevyřešených matematických problémů, které je možné formulovat pomocí naprosto elementárních prostředků. Uvažme zobrazení Collatzova hypotéza pak zní: Pro každé přirozené existuje přirozené takové, že .
nworkers()
Definujme metodu collatz
, která pro dané provede Collatzovu iteraci a vrátí počet iterací nutný k dosažení
@everywhere function collatz(n)
i = 0
while n != 1
if iseven(n)
n = div(n, 2)
else
n = 3n + 1
end
i += 1
end
return i
end
Prvních pár hodnot:
[ collatz(n) for n = 1:10 ]
Uvažme opět dvě možné paralelizace (pomocí pmap
a pomocí @distributed for
):
using SharedArrays
function collatz_pmap(nmax)
pmap(collatz, 1:nmax)
end
function collatz_dfor(nmax)
data = SharedArray{Int64}(nmax)
@sync @distributed for n = 1:nmax
data[n] = collatz(n)
end
return data
end
@time collatz_pmap(1_000_000)
@time collatz_dfor(1_000_000)
Evidentně jsou jednotlivé iterace stále dost levné, takže druhá varianta vítězí. Napočtěme si pořádnou dávku dat:
@time iterations = collatz_dfor(10^8)
using DataFrames, PyPlot
df = DataFrame(:n => 1:10^8, :iterations => iterations)
Jak vypadá distribuce počtu iterací na prvních pár číslech?
plt.xlabel(L"n")
plt.ylabel("iterations")
scatter(df.n[1:1_000_000], df.iterations[1:1_000_000], s=1)
Jaká je distribuce počtu iterací?
counts = combine(groupby(df, :iterations), nrow => :count)
plt.xlabel("iterations")
plt.ylabel("count")
scatter(counts.iterations, counts.count, s=1)
3. Profilování
Profilování pomáhá při snaze vylepšit výkon naší implementace. Oproti obyčejným benchmarkům, které jsme si ukazovali dříve, se snažíme profilováním zjistit co je podstatné pro dobu běhu programu.
V Julia lze použít modul Profile
ze standardní knihovny, který sleduje běh programu a zaznamenává kde běh "tráví nejvíce času".
Profile
periodicky zaznamenává backtrace aktuálně prováděného řádku kódu a tyto údaje agreguje.
3.1 Ukázka
Názornější bude asi jednoduchá ukázka na nějaké elementární funkci.
Uvažme následující funkci, která vyhodnocuje funkční hodnotu polynomu zadaného pomocí vektoru koeficientů [a[0], a[1], ..., a[n]]
.
Naivní implementace by byla následující:
function poly_eval_v1(coeffs::Vector{T}, x::T) where { T <: Number }
value = zero(T)
for k in 1:length(coeffs)
value += coeffs[k] * x^(k - 1)
end
return value
end
poly_eval_v1 (generic function with 1 method)
Před profilování je dobré funkci alespoň jednou vyhodnotit a tím vyvolat její kompilaci (JIT). Jinak budeme profilovat LLVM kompilátor.
poly_eval_v1([1.0, 2.0], 10.0)
21.0
Nejprve musíme importovat zmíněný modul.
using Profile
A poté k měření použít makro @profile
.
Také použijeme větší vstup, aby bylo co měřit.
Zdůrazněme, že Profile
nesleduje každý krok programu, jen po jistých časových krocích zkontroluje běh programu.
To má jistě své výhody (malé zpomalení samotného výpočtu) i nevýhody (ne zcela přesné měření).
Toto makro akumuluje výsledky i z několika běhů. Pokud chceme předchozí data odstranit, použijeme k tomu metodu Profile.clear()
.
coeffs = rand(10^8)
Profile.clear()
@profile poly_eval_v1(coeffs, -0.8945)
0.46208491495222437
K zobrazení dat o běhu programu můžeme použít hned několik nástrojů.
Základní je zobrazení přímo pomocí metody print
z modulu Profile
:
Profile.print()
Overhead ╎ [+additional indent] Count File:Line; Function ========================================================= 1╎1 @Base/math.jl:0; pow_body(x::Float64, n::Int64) ╎9047 @Base/task.jl:514; (::IJulia.var"#15#18")() ╎ 9047 ...lia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket) ╎ 9047 @Base/essentials.jl:816; invokelatest ╎ 9047 @Base/essentials.jl:819; #invokelatest#2 ╎ 9047 ...execute_request.jl:67; execute_request(socket::ZMQ.S... ╎ 9047 ...SoftGlobalScope.jl:65; softscope_include_string(m::... ╎ ╎ 9047 @Base/loading.jl:1903; include_string(mapexpr::typ... 1╎ ╎ 9047 @Base/boot.jl:370; eval 39╎ ╎ 39 In[1]:0; poly_eval_v1(coeffs::Vect... ╎ ╎ 8886 In[1]:4; poly_eval_v1(coeffs::Vect... 48╎ ╎ 48 ...e/essentials.jl:13; getindex 121╎ ╎ 121 @Base/float.jl:410; * ╎ ╎ 8717 @Base/math.jl:1165; ^ 182╎ ╎ 182 @Base/math.jl:0; pow_body(x::Float64, n::... ╎ ╎ 23 @Base/math.jl:1171; pow_body(x::Float64, n::... 23╎ ╎ ╎ 23 @Base/promotion.jl:499; == 277╎ ╎ 588 @Base/math.jl:1179; pow_body(x::Float64, n::... ╎ ╎ ╎ 311 @Base/operators.jl:369; > 311╎ ╎ ╎ 311 @Base/int.jl:83; < 55╎ ╎ 168 @Base/math.jl:1180; pow_body(x::Float64, n::... ╎ ╎ ╎ 113 @Base/operators.jl:369; > 113╎ ╎ ╎ 113 @Base/int.jl:83; < ╎ ╎ 1205 @Base/math.jl:1181; pow_body(x::Float64, n::... 163╎ ╎ ╎ 163 @Base/float.jl:410; * 1041╎ ╎ ╎ 1042 @Base/float.jl:413; muladd ╎ ╎ 1482 @Base/math.jl:1182; pow_body(x::Float64, n::... ╎ ╎ ╎ 1252 @Base/math.jl:47; two_mul 1252╎ ╎ ╎ 1252 @Base/float.jl:410; * ╎ ╎ ╎ 230 @Base/math.jl:48; two_mul ╎ ╎ ╎ 230 ...e/floatfuncs.jl:439; fma 230╎ ╎ ╎ 230 .../floatfuncs.jl:434; fma_llvm ╎ ╎ 1367 @Base/math.jl:1183; pow_body(x::Float64, n::... 1367╎ ╎ ╎ 1367 @Base/float.jl:408; + ╎ ╎ 1054 @Base/math.jl:1185; pow_body(x::Float64, n::... ╎ ╎ ╎ 1054 @Base/operators.jl:578; * 860╎ ╎ ╎ 860 @Base/float.jl:410; * ╎ ╎ ╎ 194 @Base/promotion.jl:411; * 194╎ ╎ ╎ 194 @Base/float.jl:410; * ╎ ╎ 991 @Base/math.jl:1186; pow_body(x::Float64, n::... ╎ ╎ ╎ 895 @Base/math.jl:47; two_mul 895╎ ╎ ╎ 895 @Base/float.jl:410; * ╎ ╎ ╎ 96 @Base/math.jl:48; two_mul ╎ ╎ ╎ 96 ...e/floatfuncs.jl:439; fma 96╎ ╎ ╎ 96 .../floatfuncs.jl:434; fma_llvm ╎ ╎ 104 @Base/math.jl:1187; pow_body(x::Float64, n::... 104╎ ╎ ╎ 104 @Base/float.jl:408; + ╎ ╎ 1111 @Base/math.jl:1188; pow_body(x::Float64, n::... ╎ ╎ ╎ 1111 @Base/int.jl:512; >>> 1111╎ ╎ ╎ 1111 @Base/int.jl:504; >>> ╎ ╎ 14 @Base/math.jl:1190; pow_body(x::Float64, n::... 1╎ ╎ ╎ 1 @Base/float.jl:410; * 13╎ ╎ ╎ 13 @Base/float.jl:413; muladd ╎ ╎ 428 @Base/math.jl:1191; pow_body(x::Float64, n::... 214╎ ╎ ╎ 214 ...e/essentials.jl:575; ifelse ╎ ╎ ╎ 99 @Base/float.jl:622; isfinite 99╎ ╎ ╎ 99 @Base/float.jl:409; - 115╎ ╎ ╎ 115 @Base/float.jl:413; muladd ╎ ╎ 104 In[1]:5; poly_eval_v1(coeffs::Vect... ╎ ╎ 104 @Base/range.jl:891; iterate 104╎ ╎ 104 @Base/promotion.jl:499; == 17╎ ╎ 17 @Base/math.jl:1168; pow_body(x::Float64, n::I... 1╎1 ...ia/src/heartbeat.jl:22; heartbeat_thread(sock::Ptr{Noth... Total snapshots: 9050. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.
Alternativně:
Profile.print(format=:flat)
Count Overhead File Line Function ===== ======== ==== ==== ======== 39 39 In[1] ? poly_eval_v1(coeffs::Vector{Flo... 8886 0 In[1] 4 poly_eval_v1(coeffs::Vector{Flo... 104 0 In[1] 5 poly_eval_v1(coeffs::Vector{Flo... 9047 1 @Base/boot.jl 370 eval 9047 0 @Base/essentials.jl 819 #invokelatest#2 48 48 @Base/essentials.jl 13 getindex 214 214 @Base/essentials.jl 575 ifelse 9047 0 @Base/essentials.jl 816 invokelatest 3486 3486 @Base/float.jl 410 * 1471 1471 @Base/float.jl 408 + 99 99 @Base/float.jl 409 - 99 0 @Base/float.jl 622 isfinite 1170 1169 @Base/float.jl 413 muladd 326 0 @Base/floatfuncs.jl 439 fma 326 326 @Base/floatfuncs.jl 434 fma_llvm 424 424 @Base/int.jl 83 < 1111 1111 @Base/int.jl 504 >>> 1111 0 @Base/int.jl 512 >>> 9047 0 @Base/loading.jl 1903 include_string(mapexpr::typeof(... 8717 0 @Base/math.jl 1165 ^ 183 183 @Base/math.jl ? pow_body(x::Float64, n::Int64) 17 17 @Base/math.jl 1168 pow_body(x::Float64, n::Int64) 23 0 @Base/math.jl 1171 pow_body(x::Float64, n::Int64) 588 277 @Base/math.jl 1179 pow_body(x::Float64, n::Int64) 168 55 @Base/math.jl 1180 pow_body(x::Float64, n::Int64) 1205 0 @Base/math.jl 1181 pow_body(x::Float64, n::Int64) 1482 0 @Base/math.jl 1182 pow_body(x::Float64, n::Int64) 1367 0 @Base/math.jl 1183 pow_body(x::Float64, n::Int64) 1054 0 @Base/math.jl 1185 pow_body(x::Float64, n::Int64) 991 0 @Base/math.jl 1186 pow_body(x::Float64, n::Int64) 104 0 @Base/math.jl 1187 pow_body(x::Float64, n::Int64) 1111 0 @Base/math.jl 1188 pow_body(x::Float64, n::Int64) 14 0 @Base/math.jl 1190 pow_body(x::Float64, n::Int64) 428 0 @Base/math.jl 1191 pow_body(x::Float64, n::Int64) 2147 0 @Base/math.jl 47 two_mul 326 0 @Base/math.jl 48 two_mul 1054 0 @Base/operators.jl 578 * 424 0 @Base/operators.jl 369 > 194 0 @Base/promotion.jl 411 * 127 127 @Base/promotion.jl 499 == 104 0 @Base/range.jl 891 iterate 9047 0 @Base/task.jl 514 (::IJulia.var"#15#18")() 9047 0 ...ia/src/eventloop.jl 8 eventloop(socket::ZMQ.Socket) 9047 0 .../execute_request.jl 67 execute_request(socket::ZMQ.Soc... 1 1 ...ia/src/heartbeat.jl 22 heartbeat_thread(sock::Ptr{Noth... 9047 0 .../SoftGlobalScope.jl 65 softscope_include_string(m::Mod... Total snapshots: 9050. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.
Na cvičení kontaktně podrobně rozebereme tento výpis. Zde textově by to bylo poměrně komplikované.
Úvodní část souvisí s IJulia komunikací s notebookem. Co jistě v obou výrazech vidíme, je vyčnívající řádek
5237 5234 @Base/math.jl 920 ^
Tj. náš kód tráví hodně času umocňováním.
Vedle textového znázornění můžeme použít i grafické, například pomocí balíčku ProfileSVG.jl
(tj. ] add ProfileSVG
).
using ProfileSVG
[ Info: Precompiling ProfileSVG [132c30aa-f267-4189-9183-c8a63c7e05e6]
@profview poly_eval_v1(coeffs, -0.8945)
Opět vidíme, že umocňování je nejnáročnější.
3.2 Vylepšení I
Mohli bychom tedy nějak počet umocňování zmenšit? Zcela jistě! Můžeme se ho úplně zbavit!
function poly_eval_v2(coeffs::Vector{T}, x::T) where { T <: Number }
value = zero(T)
powx = one(T)
for k in 1:length(coeffs)
value += coeffs[k] * powx
powx *= x
end
return value
end
poly_eval_v2 (generic function with 1 method)
poly_eval_v2([1., 2., 3.], 10.)
321.0
Profile.clear()
@profile poly_eval_v2(coeffs, -0.8945)
0.4620849149522243
Profile.print()
Overhead ╎ [+additional indent] Count File:Line; Function ========================================================= ╎12967 @Base/task.jl:514; (::IJulia.var"#15#18")() ╎ 12967 ...ia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket) ╎ 12967 @Base/essentials.jl:816; invokelatest ╎ 12967 @Base/essentials.jl:819; #invokelatest#2 ╎ 12967 ...execute_request.jl:67; execute_request(socket::ZMQ.... ╎ 12967 ...oftGlobalScope.jl:65; softscope_include_string(m::... ╎ ╎ 12967 @Base/loading.jl:1903; include_string(mapexpr::ty... ╎ ╎ 12967 @Base/boot.jl:370; eval ╎ ╎ 6245 In[11]:6; poly_eval_v2(coeffs::Vector... 94╎ ╎ 94 @Base/float.jl:410; * 6151╎ ╎ 6151 @Base/float.jl:408; + ╎ ╎ 602 In[11]:7; poly_eval_v2(coeffs::Vector... 602╎ ╎ 602 @Base/float.jl:410; * 6119╎ ╎ 6120 In[11]:8; poly_eval_v2(coeffs::Vector... Total snapshots: 12967. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.
@profview poly_eval_v2(coeffs, -0.8945)
Teď vidíme, že ^
vůbec nevoláme a z algebraických operací používáme jen +
a *
.
3.3 Vylepšení II
Nemohli bychom snížit i počet algebraických operací nutných k vypočtení naší funkční hodnoty? Ano, mohli! Stačí k vyhodnocování použít tzv. Hornerovu metodu, což je pouze oficiální název pro využití distributivního zákona (jiné uzávorkování). Například
"""
Vyhodnocování polynomu pomocí Hornerovy metody.
"""
function poly_eval_v3(coeffs::Vector{T}, x::T) where { T <: Number }
value = coeffs[end]
for k in (length(coeffs)-1):-1:1
value = value * x + coeffs[k]
end
return value
end
poly_eval_v3
poly_eval_v3([1., 2., 3.], 10.)
321.0
poly_eval_v3([42.], -2.)
42.0
Profile.clear()
@profile poly_eval_v3(coeffs, -0.8945)
0.46208491495222453
Profile.print()
Overhead ╎ [+additional indent] Count File:Line; Function ========================================================= ╎267 @Base/task.jl:514; (::IJulia.var"#15#18")() ╎ 267 @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket) ╎ 267 @Base/essentials.jl:816; invokelatest ╎ 267 @Base/essentials.jl:819; #invokelatest#2 ╎ 267 .../execute_request.jl:67; execute_request(socket::ZMQ.So... ╎ 267 .../SoftGlobalScope.jl:65; softscope_include_string(m::Mo... ╎ ╎ 267 @Base/loading.jl:1903; include_string(mapexpr::type... ╎ ╎ 267 @Base/boot.jl:370; eval ╎ ╎ 84 In[16]:7; poly_eval_v3(coeffs::Vector{... 84╎ ╎ 84 @Base/float.jl:408; + ╎ ╎ 183 In[16]:8; poly_eval_v3(coeffs::Vector{... ╎ ╎ 183 @Base/range.jl:891; iterate 182╎ ╎ 183 @Base/promotion.jl:499; == Total snapshots: 267. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.
@profview poly_eval_v3(coeffs, -0.8945)
Vidíme, že násobení i sčítání proporcionálně jistě zabírá méně času než v předchozím případě.
Pokud se podíváme na profil tak nás asi zaujmou poslední dva řádky:
╎ ╎ 258 @Base/range.jl:674; iterate
257╎ ╎ 258 @Base/promotion.jl:410; ==
Tyto má dle profilu na svědomí řádek 7 v buňce kde definujeme poly_eval_v3
, tj.
value = value * x + coeffs[k]
Co se zde děje? Podle popisu půjde spíše o efekt související s iterátorem, resp. for
cyklem.
Nemohli bychom iterovat efektivněji?
"""
Vyhodnocování polynomu pomocí Hornerovy metody, alternativní iterování.
"""
function poly_eval_v4(reversed_coeffs::Vector{T}, x::T) where { T <: Number }
value = reversed_coeffs[1]
for k = 2:length(reversed_coeffs)
value = value * x + reversed_coeffs[k]
end
return value
end
poly_eval_v4
poly_eval_v4([3., 2., 1.], 10.)
321.0
Profile.clear()
reversed_coeffs = reverse(coeffs)
@profile poly_eval_v4(reversed_coeffs, -0.8945)
0.46208491495222453
Profile.print()
Overhead ╎ [+additional indent] Count File:Line; Function ========================================================= ╎253 @Base/task.jl:514; (::IJulia.var"#15#18")() ╎ 253 @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket) ╎ 253 @Base/essentials.jl:816; invokelatest ╎ 253 @Base/essentials.jl:819; #invokelatest#2 ╎ 253 .../execute_request.jl:67; execute_request(socket::ZMQ.So... ╎ 253 .../SoftGlobalScope.jl:65; softscope_include_string(m::Mo... ╎ ╎ 253 @Base/loading.jl:1903; include_string(mapexpr::type... ╎ ╎ 253 @Base/boot.jl:370; eval ╎ ╎ 104 In[22]:7; poly_eval_v4(reversed_coeffs... 103╎ ╎ 104 @Base/float.jl:408; + 149╎ ╎ 149 In[22]:8; poly_eval_v4(reversed_coeffs... Total snapshots: 253. Utilization: 100% across all threads and tasks. Use the `groupby` kwarg to break down by thread and/or task.
@profview poly_eval_v4(coeffs, -0.8945)
Zdá se, že jsme lehce ubrali. Zkusme jednotlivé implementace nyní porovnat.
using BenchmarkTools
coeffs = rand(10^6);
@benchmark poly_eval_v1($coeffs, -0.5)
BenchmarkTools.Trial: 33 samples with 1 evaluation. Range (min … max): 151.220 ms … 167.929 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 153.338 ms ┊ GC (median): 0.00% Time (mean ± σ): 154.340 ms ± 3.044 ms ┊ GC (mean ± σ): 0.00% ± 0.00% █▁▄▁█ ▁▁▁▁ ▁ ▆▁▆█████▁████▆▁▁▆█▁▆▁▁▁▁▁▆▁▁▁▁▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁ 151 ms Histogram: frequency by time 168 ms < Memory estimate: 0 bytes, allocs estimate: 0.
@benchmark poly_eval_v2($coeffs, -0.5)
BenchmarkTools.Trial: 3004 samples with 1 evaluation. Range (min … max): 1.404 ms … 3.080 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 1.573 ms ┊ GC (median): 0.00% Time (mean ± σ): 1.650 ms ± 234.732 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ██▃▁ ▄████▇██▇▆▅▆▆▅▅▅▄▄▄▄▄▃▃▃▄▃▃▂▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁ ▃ 1.4 ms Histogram: frequency by time 2.43 ms < Memory estimate: 0 bytes, allocs estimate: 0.
@benchmark poly_eval_v3($coeffs, -0.5)
BenchmarkTools.Trial: 1876 samples with 1 evaluation. Range (min … max): 2.409 ms … 3.993 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 2.573 ms ┊ GC (median): 0.00% Time (mean ± σ): 2.651 ms ± 241.257 μs ┊ GC (mean ± σ): 0.00% ± 0.00% █▂▂▂▁▁ ▁ ████████▇██▇▅▆▆▅▅▄▄▃▄▄▃▄▄▄▃▃▃▃▃▃▃▃▃▃▃▃▂▃▃▃▃▂▂▃▃▂▂▂▂▂▂▂▂▂▂▁▂ ▃ 2.41 ms Histogram: frequency by time 3.46 ms < Memory estimate: 0 bytes, allocs estimate: 0.
@benchmark poly_eval_v4($coeffs, -0.5)
BenchmarkTools.Trial: 1917 samples with 1 evaluation. Range (min … max): 2.385 ms … 3.572 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 2.529 ms ┊ GC (median): 0.00% Time (mean ± σ): 2.594 ms ± 209.906 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▆█ ██▆███▅▆▅▆▆▅▄▆▄▅▅▄▃▃▃▃▃▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▂▁▂▂▂▁▁▁▂▂▁▁▁▁▁▁ ▂ 2.39 ms Histogram: frequency by time 3.25 ms < Memory estimate: 0 bytes, allocs estimate: 0.
"""
Vyhodnocování polynomu pomocí Hornerovy metody, alternativní iterování.
"""
function poly_eval_final(reversed_coeffs::Vector{T}, x::T) where { T <: Number }
value = reversed_coeffs[1]
for k = 2:length(reversed_coeffs)
value = muladd(value, x, reversed_coeffs[k])
end
return value
end
poly_eval_final
@benchmark poly_eval_final($coeffs, -0.5)
BenchmarkTools.Trial: 3550 samples with 1 evaluation. Range (min … max): 1.226 ms … 2.440 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 1.336 ms ┊ GC (median): 0.00% Time (mean ± σ): 1.396 ms ± 167.426 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▂█▅▁ ▅████▆█▇▇▆▅▅▄▄▄▄▄▄▄▄▃▃▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▂▁▂▁▂▁▁▁▁▂▁▁▁▁ ▂ 1.23 ms Histogram: frequency by time 1.94 ms < Memory estimate: 0 bytes, allocs estimate: 0.
"""
Vyhodnocování polynomu pomocí Hornerovy metody, alternativní iterování.
"""
function poly_eval_v6(reversed_coeffs::Vector{T}, x::T) where { T <: Number }
value = reversed_coeffs[1]
for k = 2:length(reversed_coeffs)
value *= x
value += reversed_coeffs[k]
end
return value
end
poly_eval_v6
@benchmark poly_eval_v6($coeffs, -0.5)
BenchmarkTools.Trial: 1849 samples with 1 evaluation. Range (min … max): 2.394 ms … 3.925 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 2.602 ms ┊ GC (median): 0.00% Time (mean ± σ): 2.687 ms ± 271.034 μs ┊ GC (mean ± σ): 0.00% ± 0.00% █▁▂▃ ▆█████▇▆█▇▇▅▆▅▅▅▅▄▅▅▃▄▃▃▃▃▃▃▃▃▂▂▂▂▂▂▃▃▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▂▁▁▂ ▃ 2.39 ms Histogram: frequency by time 3.56 ms < Memory estimate: 0 bytes, allocs estimate: 0.
3.4 Poznámka
Nástrojů pro analýzu profilu je více. Zajímavou alternativou je PProf.jl
založený na externím programu pprof
od Google.
Řešení některých úkolů
"""
Vyhodnocování polynomu bez umocňování.
"""
function poly_eval_v2(coeffs::Vector{T}, x::T) where { T <: Number }
value = zero(T)
power = one(T)
for k in 1:length(coeffs)
value += coeffs[k] * power
power *= x
end
return value
end
"""
Vyhodnocování polynomu pomocí Hornerovy metody.
"""
function poly_eval_v3(coeffs::Vector{T}, x::T) where { T <: Number }
value = coeffs[end]
for k in (length(coeffs)-1):-1:1
value = value * x + coeffs[k]
end
return value
end
@everywhere function collatz(n)
k = 0
while n != 1
if isodd(n)
n = 3n+1
else
n = div(n, 2)
end
k += 1
end
return k
end
Reference
Announcing composable multi-threaded parallelism in Julia
Profilování se věnuje samostatná stránka Julia dokumentace Profiling.
O distribuovaných výpočtech se více dočtene na stránce Multi-processing and Distributed Computing.