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 = 13
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 = 13
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 = 100
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 = 100
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
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.
"""
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.
"""
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) .- 10);
@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, 1000)
ims = LinRange(-1, 1, 600)
data = zeros(Int64, length(res), length(ims))
julia_matrix!(-0.8 + 0.156im, 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.156im, 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, 1000)
ims = LinRange(-1, 1, 600)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p1!(-0.8 + 0.156im, 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.156im, 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, 1000)
ims = LinRange(-1, 1, 600)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p2!(-0.8 + 0.156im, 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.156im, 10.0, $res, $ims, $data)
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, 1000)
ims = LinRange(-1, 1, 600)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p3!(-0.8 + 0.156im, 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.156im, 10.0, $res, $ims, $data)
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.
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
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, 1000)
Kolik času nám zabere procházek o stejném počtu kroků?
@time get_walks(10_000, 10_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.
function get_walks_p1(n, steps)
output = fill(steps, n)
output = pmap(random_walk, output)
return output
end
get_walks_p1(10, 1000)
Kolik času nám zabere procházek o stejném počtu kroků?
@time get_walks_p1(10_000, 10_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, 1000)
Kolik času nám zabere procházek o stejném počtu kroků?
@time get_walks_p2(10_000, 10_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) # "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 + 1000)
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()
.
Profile.clear()
coeffs = rand(10^8)
@profile poly_eval_v1(coeffs, -0.8945)
0.3855995607565534
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 ========================================================= ╎9213 @Base/task.jl:484; (::IJulia.var"#15#18")() ╎ 9213 ...lia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket) ╎ 9213 @Base/essentials.jl:726; invokelatest ╎ 9213 @Base/essentials.jl:729; #invokelatest#2 ╎ 9213 ...execute_request.jl:67; execute_request(socket::ZMQ.S... ╎ 9213 ...SoftGlobalScope.jl:65; softscope_include_string(m::... ╎ ╎ 9213 @Base/loading.jl:1428; include_string(mapexpr::typ... 1╎ ╎ 9213 @Base/boot.jl:368; eval ╎ ╎ 9086 In[1]:4; poly_eval_v1(coeffs::Vect... 37╎ ╎ 37 @Base/array.jl:924; getindex 27╎ ╎ 27 @Base/float.jl:385; * 79╎ ╎ 79 @Base/float.jl:383; + 23╎ ╎ 23 @Base/math.jl:0; ^ 32╎ ╎ 8920 @Base/math.jl:1046; ^ 122╎ ╎ 122 @Base/math.jl:0; pow_body(x::Float64, n::... ╎ ╎ 1 @Base/math.jl:1053; pow_body(x::Float64, n::... 1╎ ╎ ╎ 1 @Base/int.jl:83; < 315╎ ╎ 600 @Base/math.jl:1060; pow_body(x::Float64, n::... ╎ ╎ ╎ 285 @Base/operators.jl:382; > 285╎ ╎ ╎ 285 @Base/int.jl:83; < ╎ ╎ 305 @Base/math.jl:1061; pow_body(x::Float64, n::... ╎ ╎ ╎ 305 @Base/operators.jl:382; > 305╎ ╎ ╎ 305 @Base/int.jl:83; < ╎ ╎ 1482 @Base/math.jl:1062; pow_body(x::Float64, n::... 177╎ ╎ ╎ 177 @Base/float.jl:385; * 1305╎ ╎ ╎ 1305 @Base/float.jl:388; muladd ╎ ╎ 1393 @Base/math.jl:1063; pow_body(x::Float64, n::... ╎ ╎ ╎ 1155 @Base/math.jl:47; two_mul 1155╎ ╎ ╎ 1155 @Base/float.jl:385; * ╎ ╎ ╎ 238 @Base/math.jl:48; two_mul ╎ ╎ ╎ 238 ...e/floatfuncs.jl:422; fma 238╎ ╎ ╎ 238 .../floatfuncs.jl:417; fma_llvm ╎ ╎ 1253 @Base/math.jl:1064; pow_body(x::Float64, n::... 1253╎ ╎ ╎ 1253 @Base/float.jl:383; + ╎ ╎ 984 @Base/math.jl:1066; pow_body(x::Float64, n::... ╎ ╎ ╎ 984 @Base/operators.jl:591; * 793╎ ╎ ╎ 793 @Base/float.jl:385; * ╎ ╎ ╎ 191 @Base/promotion.jl:389; * 191╎ ╎ ╎ 191 @Base/float.jl:385; * ╎ ╎ 1204 @Base/math.jl:1067; pow_body(x::Float64, n::... ╎ ╎ ╎ 1023 @Base/math.jl:47; two_mul 1023╎ ╎ ╎ 1023 @Base/float.jl:385; * ╎ ╎ ╎ 181 @Base/math.jl:48; two_mul ╎ ╎ ╎ 181 ...e/floatfuncs.jl:422; fma 181╎ ╎ ╎ 181 .../floatfuncs.jl:417; fma_llvm ╎ ╎ 76 @Base/math.jl:1068; pow_body(x::Float64, n::... 76╎ ╎ ╎ 76 @Base/float.jl:383; + ╎ ╎ 1124 @Base/math.jl:1069; pow_body(x::Float64, n::... ╎ ╎ ╎ 1124 @Base/int.jl:505; >>> 1123╎ ╎ ╎ 1124 @Base/int.jl:497; >>> ╎ ╎ 7 @Base/math.jl:1071; pow_body(x::Float64, n::... 7╎ ╎ ╎ 7 @Base/float.jl:388; muladd ╎ ╎ 337 @Base/math.jl:1072; pow_body(x::Float64, n::... 84╎ ╎ ╎ 84 ...e/essentials.jl:489; ifelse ╎ ╎ ╎ 148 @Base/float.jl:499; isfinite 148╎ ╎ ╎ 148 @Base/float.jl:384; - 105╎ ╎ ╎ 105 @Base/float.jl:388; muladd ╎ ╎ 124 In[1]:5; poly_eval_v1(coeffs::Vect... ╎ ╎ 124 @Base/range.jl:883; iterate 124╎ ╎ 124 @Base/promotion.jl:477; == ╎ ╎ 2 @Base/math.jl:1052; pow_body(x::Float64, n::I... 2╎ ╎ 2 @Base/promotion.jl:477; == Total snapshots: 9215. 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 ===== ======== ==== ==== ======== 9086 0 In[1] 4 poly_eval_v1(coeffs::Vector{Flo... 124 0 In[1] 5 poly_eval_v1(coeffs::Vector{Flo... 37 0 @Base/array.jl 924 getindex 9213 0 @Base/boot.jl 368 eval 9213 0 @Base/essentials.jl 729 #invokelatest#2 84 0 @Base/essentials.jl 489 ifelse 9213 0 @Base/essentials.jl 726 invokelatest 3366 0 @Base/float.jl 385 * 1408 0 @Base/float.jl 383 + 148 0 @Base/float.jl 384 - 148 0 @Base/float.jl 499 isfinite 1417 0 @Base/float.jl 388 muladd 419 0 @Base/floatfuncs.jl 422 fma 419 0 @Base/floatfuncs.jl 417 fma_llvm 591 0 @Base/int.jl 83 < 1124 0 @Base/int.jl 497 >>> 1124 0 @Base/int.jl 505 >>> 9213 0 @Base/loading.jl 1428 include_string(mapexpr::typeof(... 23 0 @Base/math.jl ? ^ 8920 0 @Base/math.jl 1046 ^ 122 0 @Base/math.jl ? pow_body(x::Float64, n::Int64) 2 0 @Base/math.jl 1052 pow_body(x::Float64, n::Int64) 1 0 @Base/math.jl 1053 pow_body(x::Float64, n::Int64) 600 0 @Base/math.jl 1060 pow_body(x::Float64, n::Int64) 305 0 @Base/math.jl 1061 pow_body(x::Float64, n::Int64) 1482 0 @Base/math.jl 1062 pow_body(x::Float64, n::Int64) 1393 0 @Base/math.jl 1063 pow_body(x::Float64, n::Int64) 1253 0 @Base/math.jl 1064 pow_body(x::Float64, n::Int64) 984 0 @Base/math.jl 1066 pow_body(x::Float64, n::Int64) 1204 0 @Base/math.jl 1067 pow_body(x::Float64, n::Int64) 76 0 @Base/math.jl 1068 pow_body(x::Float64, n::Int64) 1124 0 @Base/math.jl 1069 pow_body(x::Float64, n::Int64) 7 0 @Base/math.jl 1071 pow_body(x::Float64, n::Int64) 337 0 @Base/math.jl 1072 pow_body(x::Float64, n::Int64) 2178 0 @Base/math.jl 47 two_mul 419 0 @Base/math.jl 48 two_mul 984 0 @Base/operators.jl 591 * 590 0 @Base/operators.jl 382 > 191 0 @Base/promotion.jl 389 * 126 0 @Base/promotion.jl 477 == 124 0 @Base/range.jl 883 iterate 9213 9213 @Base/task.jl 484 (::IJulia.var"#15#18")() 9213 0 ...ia/src/eventloop.jl 8 eventloop(socket::ZMQ.Socket) 9213 0 .../execute_request.jl 67 execute_request(socket::ZMQ.Soc... 9213 0 .../SoftGlobalScope.jl 65 softscope_include_string(m::Mod... Total snapshots: 9215 (100% utilization 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
@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.3855995607565533
Profile.print()
Overhead ╎ [+additional indent] Count File:Line; Function ========================================================= ╎7852 @Base/task.jl:484; (::IJulia.var"#15#18")() ╎ 7852 ...lia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket) ╎ 7852 @Base/essentials.jl:726; invokelatest ╎ 7852 @Base/essentials.jl:729; #invokelatest#2 ╎ 7852 ...execute_request.jl:67; execute_request(socket::ZMQ.S... ╎ 7852 ...SoftGlobalScope.jl:65; softscope_include_string(m::... ╎ ╎ 7852 @Base/loading.jl:1428; include_string(mapexpr::typ... ╎ ╎ 7852 @Base/boot.jl:368; eval ╎ ╎ 3857 In[9]:6; poly_eval_v2(coeffs::Vector... 63╎ ╎ 63 @Base/float.jl:385; * 3794╎ ╎ 3794 @Base/float.jl:383; + ╎ ╎ 386 In[9]:7; poly_eval_v2(coeffs::Vector... 386╎ ╎ 386 @Base/float.jl:385; * 3608╎ ╎ 3609 In[9]:8; poly_eval_v2(coeffs::Vector... Total snapshots: 7852. 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.3855995607565532
Profile.print()
Overhead ╎ [+additional indent] Count File:Line; Function ========================================================= ╎238 @Base/task.jl:484; (::IJulia.var"#15#18")() ╎ 238 @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket) ╎ 238 @Base/essentials.jl:726; invokelatest ╎ 238 @Base/essentials.jl:729; #invokelatest#2 ╎ 238 .../execute_request.jl:67; execute_request(socket::ZMQ.So... ╎ 238 .../SoftGlobalScope.jl:65; softscope_include_string(m::Mo... ╎ ╎ 238 @Base/loading.jl:1428; include_string(mapexpr::type... ╎ ╎ 238 @Base/boot.jl:368; eval ╎ ╎ 72 In[14]:7; poly_eval_v3(coeffs::Vector{... 72╎ ╎ 72 @Base/float.jl:383; + ╎ ╎ 166 In[14]:8; poly_eval_v3(coeffs::Vector{... ╎ ╎ 166 @Base/range.jl:883; iterate 165╎ ╎ 166 @Base/promotion.jl:477; == Total snapshots: 238. 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.3855995607565532
Profile.print()
Overhead ╎ [+additional indent] Count File:Line; Function ========================================================= ╎224 @Base/task.jl:484; (::IJulia.var"#15#18")() ╎ 224 @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket) ╎ 224 @Base/essentials.jl:726; invokelatest ╎ 224 @Base/essentials.jl:729; #invokelatest#2 ╎ 224 .../execute_request.jl:67; execute_request(socket::ZMQ.So... ╎ 224 .../SoftGlobalScope.jl:65; softscope_include_string(m::Mo... ╎ ╎ 224 @Base/loading.jl:1428; include_string(mapexpr::type... ╎ ╎ 224 @Base/boot.jl:368; eval ╎ ╎ 109 In[20]:7; poly_eval_v4(reversed_coeffs... 108╎ ╎ 109 @Base/float.jl:383; + 115╎ ╎ 115 In[20]:8; poly_eval_v4(reversed_coeffs... Total snapshots: 224. 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.940 ms … 153.603 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 152.334 ms ┊ GC (median): 0.00% Time (mean ± σ): 152.385 ms ± 319.044 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▁▁ ▃ █ ▄▁▄▄▁▄▄██▄▁▄▄▄▇▁▁▇▄█▁▁█▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▁ 152 ms Histogram: frequency by time 154 ms < Memory estimate: 0 bytes, allocs estimate: 0.
@benchmark poly_eval_v2($coeffs, -0.5)
BenchmarkTools.Trial: 3565 samples with 1 evaluation. Range (min … max): 1.392 ms … 1.534 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 1.398 ms ┊ GC (median): 0.00% Time (mean ± σ): 1.399 ms ± 8.413 μs ┊ GC (mean ± σ): 0.00% ± 0.00% █▆ ▂▆▃ ▆███▇█████▇▅▄▄▄▄▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▂▁▂▂▂▂▁▂▂▂▂ ▃ 1.39 ms Histogram: frequency by time 1.44 ms < Memory estimate: 0 bytes, allocs estimate: 0.
@benchmark poly_eval_v3($coeffs, -0.5)
BenchmarkTools.Trial: 2083 samples with 1 evaluation. Range (min … max): 2.384 ms … 2.536 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 2.396 ms ┊ GC (median): 0.00% Time (mean ± σ): 2.397 ms ± 9.233 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▁▂▆███▅▇▅▃▄▃▁ ▁ ▂▂▃▄▇▇▆████████████████▇▇▆▅▅▅▅▄▄▃▃▂▃▃▃▂▃▂▃▃▂▂▂▂▂▂▂▂▂▁▁▁▁▂ ▄ 2.38 ms Histogram: frequency by time 2.43 ms < Memory estimate: 0 bytes, allocs estimate: 0.
@benchmark poly_eval_v4($coeffs, -0.5)
BenchmarkTools.Trial: 2101 samples with 1 evaluation. Range (min … max): 2.368 ms … 2.462 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 2.376 ms ┊ GC (median): 0.00% Time (mean ± σ): 2.377 ms ± 5.406 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▄█▇█▁ ▁▄▇▇█▅▁ ▂▂▃▅▇▆▅▄▅▅█████████████▇▇▅▅▅▄▄▃▃▃▂▂▂▃▃▂▂▂▂▁▂▂▂▂▁▂▁▂▂▁▁▁▁▂ ▄ 2.37 ms Histogram: frequency by time 2.4 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: 4074 samples with 1 evaluation. Range (min … max): 1.216 ms … 1.307 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 1.223 ms ┊ GC (median): 0.00% Time (mean ± σ): 1.224 ms ± 6.370 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▄█▇▅▂▁▁▂▃▂ ▂▃▃████████████▇▆▅▄▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂ ▃ 1.22 ms Histogram: frequency by time 1.26 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: 2102 samples with 1 evaluation. Range (min … max): 2.368 ms … 2.448 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 2.375 ms ┊ GC (median): 0.00% Time (mean ± σ): 2.376 ms ± 6.059 μs ┊ GC (mean ± σ): 0.00% ± 0.00% ▄▅▄▂▄█▅▂ ▂▄▅▄▇████████▆▅▄▃▃▃▃▃▃▂▃▂▂▂▂▂▂▁▁▂▂▂▂▂▂▂▁▂▂▁▂▂▁▂▂▁▂▁▂▂▂▂▁▂ ▃ 2.37 ms Histogram: frequency by time 2.41 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.