Jdi na navigaci předmětu

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()

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, skylake)
Threads: 4 default, 0 interactive, 2 GC (on 8 virtual cores)
Environment:
  JULIA_NUM_THREADS = 4

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

10-element Vector{Int64}:
 1
 1
 1
 1
 1
 1
 1
 1
 1
 1

import .Threads: @threads

# pole o deseti prvcích
L   = 10
zrs = zeros(Int64, L)

@threads :dynamic for j in 1:L
    # do každého prvku zapíšeme id vlákna
    zrs[j] = Threads.threadid()
end

zrs

10-element Vector{Int64}:
 1
 1
 1
 2
 2
 2
 3
 3
 4
 4

# pole o deseti prvcích
L   = 10
zrs = zeros(Int64, L)

@threads :static for j in 1:L
    # do každého prvku zapíšeme id vlákna
    zrs[j] = Threads.threadid()
end

zrs

10-element Vector{Int64}:
 1
 1
 1
 2
 2
 2
 3
 3
 4
 4

# pole o deseti prvcích
L   = 10
zrs = zeros(Int64, L)

@threads :greedy for j in 1:L
    # do každého prvku zapíšeme id vlákna
    zrs[j] = Threads.threadid()
end

zrs

10-element Vector{Int64}:
 4
 2
 2
 3
 2
 2
 3
 1
 2
 3

Co když každý krok trvá jinou dobu?

# pole o deseti prvcích
L   = 20
zrs = zeros(Int64, L)

@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

20-element Vector{Int64}:
 3
 3
 3
 3
 3
 2
 2
 2
 1
 1
 3
 3
 3
 1
 2
 4
 4
 4
 4
 4

# pole o deseti prvcích
L   = 20
zrs = zeros(Int64, L)

@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

20-element Vector{Int64}:
 1
 1
 1
 1
 1
 2
 2
 2
 2
 2
 3
 3
 3
 3
 3
 4
 4
 4
 4
 4

# pole o deseti prvcích
L   = 20
zrs = zeros(Int64, L)

@threads :greedy for j in 1:L
    sleep(rand()/2)
    # do každého prvku zapíšeme id vlákna
    zrs[j] = Threads.threadid()
end

zrs

20-element Vector{Int64}:
 2
 3
 1
 2
 1
 3
 3
 2
 2
 2
 2
 3
 4
 4
 4
 3
 1
 4
 2
 1

`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

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

f (generic function with 1 method)

Jaký je výsledek?

f()

Kolik času to zabralo?

@benchmark f()

BenchmarkTools.Trial: 1 sample with 1 evaluation.
 Single result which took 8.004 s (0.00% GC) to evaluate,
 with a memory estimate of 5.00 KiB, over 185 allocations.

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 :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 :static for j = 1:Threads.nthreads()
        values[j] = sleepy_function()
    end

    return sum(values)
end

"""
@threads verze (:greedy).
"""
function i()
    values = zeros(Int64, Threads.nthreads())

    @threads :greedy for j = 1:Threads.nthreads()
        values[j] = sleepy_function()
    end

    return sum(values)
end;

g()

h()

ch()

i()

@benchmark f()

BenchmarkTools.Trial: 1 sample with 1 evaluation.
 Single result which took 8.005 s (0.00% GC) to evaluate,
 with a memory estimate of 5.00 KiB, over 185 allocations.

@benchmark g()

BenchmarkTools.Trial: 3 samples with 1 evaluation.
 Range (min … max):  2.001 s …    2.001 s  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.001 s               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.001 s ± 450.027 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █                        █                               █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  2 s            Histogram: frequency by time            2 s <

 Memory estimate: 3.95 KiB, allocs estimate: 85.

@benchmark h()

BenchmarkTools.Trial: 3 samples with 1 evaluation.
 Range (min … max):  2.001 s …    2.001 s  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.001 s               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.001 s ± 147.355 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █         █                                              █  
  █▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  2 s            Histogram: frequency by time            2 s <

 Memory estimate: 3.59 KiB, allocs estimate: 81.

@benchmark ch()

BenchmarkTools.Trial: 3 samples with 1 evaluation.
 Range (min … max):  2.001 s …    2.002 s  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.001 s               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.001 s ± 440.552 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █                      █                                 █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  2 s            Histogram: frequency by time            2 s <

 Memory estimate: 3.59 KiB, allocs estimate: 81.

@benchmark i()

BenchmarkTools.Trial: 3 samples with 1 evaluation.
 Range (min … max):  2.001 s …    2.002 s  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.001 s               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.001 s ± 561.796 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █                   █                                    █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  2 s            Histogram: frequency by time            2 s <

 Memory estimate: 4.88 KiB, allocs estimate: 100.

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 :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 :static for j in axes(m, 2)
        @inbounds m[target, j] += α * m[source, j]
    end
end

function add9!(m::Matrix{T}, α::T, source::Integer, target::Integer) where { T <: Number }
    @threads :greedy for j in axes(m, 2)
        @inbounds m[target, j] += α * m[source, j]
    end
end;

using BenchmarkTools

mat = 20 * (rand(4, 10^6) .- 0.5);

@benchmark add1!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 273 samples with 1 evaluation.
 Range (min … max):   9.727 ms … 77.584 ms  ┊ GC (min … max): 4.98% … 79.17%
 Time  (median):     13.162 ms              ┊ GC (median):    4.15%
 Time  (mean ± σ):   13.640 ms ±  4.477 ms  ┊ GC (mean ± σ):  5.93% ±  4.83%

         ▂ ▁ ▃▅▅▅▇▂▇▄█▅▁▁                                      
  ▅▃▆▄▆▆▇█▆█▆████████████▅▃▄▃▄▃▆▃▄▃▃▁▁▃▁▃▃▃▁▃▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▃ ▃
  9.73 ms         Histogram: frequency by time        22.6 ms <

 Memory estimate: 30.52 MiB, allocs estimate: 12.

@benchmark add2!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 273 samples with 1 evaluation.
 Range (min … max):   9.178 ms … 72.932 ms  ┊ GC (min … max): 0.00% … 81.99%
 Time  (median):     13.391 ms              ┊ GC (median):    4.07%
 Time  (mean ± σ):   13.774 ms ±  5.426 ms  ┊ GC (mean ± σ):  7.36% ±  7.01%

                     ▃ ▁▇▆▄▅█  ▁                               
  ▃▅▄▇▄▅▃▃▅▆▆▆▆▆▆▆▄███████████▇█▃▆▃▃▃▅▃▅▆▃▄▃▃▃▁▃▄▃▃▁▃▃▁▃▁▁▁▁▃ ▃
  9.18 ms         Histogram: frequency by time          20 ms <

 Memory estimate: 30.52 MiB, allocs estimate: 12.

@benchmark add3!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 472 samples with 1 evaluation.
 Range (min … max):  2.342 ms … 8.742 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.335 ms             ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.579 ms ± 1.033 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅█▃▃  ▁                                                    
  ███████▇▆▅▆▄█▆▄▄▅▄▅▄▃▄▂▃▄▇▇▆▅▄▃▆▅▄▃▅▅▂▂▄▄▃▂▂▂▃▃▃▄▃▂▄▁▃▁▂▂ ▃
  2.34 ms        Histogram: frequency by time       5.87 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark add4!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 507 samples with 1 evaluation.
 Range (min … max):  2.364 ms …   6.782 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.999 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.318 ms ± 856.614 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▄   ▁▁ ▁                                                    
  ███▆▄████▆▅▅▄▄▅▃▃▆▄▆▆▄▃▃▅▃▄▄▄▃▅▄▃▄▃▂▃▄▃▄▄▄▄▅▄▃▃▁▃▃▃▂▂▃▃▂▂▁▂ ▃
  2.36 ms         Histogram: frequency by time        5.46 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark add5!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 504 samples with 1 evaluation.
 Range (min … max):  2.357 ms …   6.212 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.086 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.321 ms ± 828.574 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▇▁  ▂▃ ▃                                                    
  ███▇██████▄▅▅▃▃▆▄▃▄▄▆▄▅▅▇▂▅▄▄▃▄▃▆▃▄▄▅▄▃▅▃▄▄▆▇█▅▃▄▃▄▃▃▃▁▂▁▂▂ ▄
  2.36 ms         Histogram: frequency by time        5.19 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark add6!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 503 samples with 1 evaluation.
 Range (min … max):  2.366 ms …   7.493 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.129 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.383 ms ± 927.161 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▆                                                           
  ███▄▅▇▅▅▆▄▄▂▂▃▃▄▃▃▃▃▄▄▄▃▂▃▃▃▂▃▄▃▃▃▂▃▃▂▂▄▄▂▂▂▂▄▂▂▃▂▂▂▂▂▁▁▂▂▁ ▂
  2.37 ms         Histogram: frequency by time        5.57 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark add7!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 667 samples with 1 evaluation.
 Range (min … max):  2.403 ms …   5.119 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     3.085 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.122 ms ± 507.349 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅█▅        ▁            ▁       ▁                           
  ▄███▇▄▆▃▅▆▇▇█▇▇▅▃▅▆▆▅▅▄▅██▆▃▇▆█▇▇█▇▅▂▃▅▅▃▂▃▃▃▂▃▂▂▂▁▂▂▂▁▁▁▁▂ ▃
  2.4 ms          Histogram: frequency by time        4.43 ms <

 Memory estimate: 2.22 KiB, allocs estimate: 22.

@benchmark add8!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 539 samples with 1 evaluation.
 Range (min … max):  2.408 ms …   5.169 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.986 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.082 ms ± 498.998 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▁█      ▃                                                   
  ▅███▅▄▃▄▆█▆▆█▅▆▆▆▅▄▄▆▄▆▅▅▄▅▆▆▆▆▅▆▅▆▅▄▃▁▃▂▄▃▂▃▂▃▃▂▃▃▃▁▁▂▂▂▁▃ ▃
  2.41 ms         Histogram: frequency by time        4.48 ms <

 Memory estimate: 2.22 KiB, allocs estimate: 22.

@benchmark add9!(a, -2.376, 1, 2) setup=(a=copy(mat))

BenchmarkTools.Trial: 2 samples with 1 evaluation.
 Range (min … max):  2.679 s …   2.695 s  ┊ GC (min … max): 0.00% … 0.04%
 Time  (median):     2.687 s              ┊ GC (median):    0.02%
 Time  (mean ± σ):   2.687 s ± 11.546 ms  ┊ GC (mean ± σ):  0.02% ± 0.03%

  █                                                       █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  2.68 s         Histogram: frequency by time         2.7 s <

 Memory estimate: 15.26 MiB, allocs estimate: 999630.

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

julia_matrix!

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]])

PyObject <matplotlib.image.AxesImage object at 0x7fe660327380>

A podívejme se na dobu běhu.

@benchmark julia_matrix!(-0.8 + 0.256im, 10.0, $res, $ims, $data)

BenchmarkTools.Trial: 28 samples with 1 evaluation.
 Range (min … max):  168.124 ms … 222.016 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     178.220 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   180.290 ms ±  10.715 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▃▃ ▃▃▃ ▃▃ ▃ █                                               
  ▇▁▁██▁███▁██▇█▇█▇▁▁▇▁▇▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ ▁
  168 ms           Histogram: frequency by time          222 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

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 $f$ ) 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 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

julia_matrix_p1! (generic function with 1 method)

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]])

PyObject <matplotlib.image.AxesImage object at 0x7fe6608ad130>

@benchmark julia_matrix_p1!(-0.8 + 0.256im, 10.0, $res, $ims, $data)

BenchmarkTools.Trial: 72 samples with 1 evaluation.
 Range (min … max):  62.205 ms … 134.551 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     65.494 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   70.335 ms ±  11.585 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅█▄                                                          
  ▇███▇█▃▃▃▃▁▁▁▃▃▁▁▃▄▇▃▄▃▃▁▃▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▁
  62.2 ms         Histogram: frequency by time          114 ms <

 Memory estimate: 2.47 KiB, allocs estimate: 22.

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 :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

julia_matrix_p2! (generic function with 1 method)

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]])

PyObject <matplotlib.image.AxesImage object at 0x7fe6603cc950>

@benchmark julia_matrix_p2!(-0.8 + 0.256im, 10.0, $res, $ims, $data)

BenchmarkTools.Trial: 60 samples with 1 evaluation.
 Range (min … max):  64.386 ms … 174.248 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     79.906 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   83.567 ms ±  18.429 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁▄█▁██      ▁  ▁ ▁▄▁          ▁  ▁   ▁  ▁   ▄                 
  ██████▁▆▆▆▁▁█▁▁█▆███▁▆▆▆▆▁▆▁▁▆█▁▁█▁▆▆█▁▆█▁▁▆█▁▁▁▆▁▁▁▁▁▆▁▆▁▁▆ ▁
  64.4 ms         Histogram: frequency by time          115 ms <

 Memory estimate: 2.47 KiB, allocs estimate: 22.

function julia_matrix_p3!(c::Complex{Float64}, R::Float64, reals, imags, data; imax::Int64=1000)
    f(u) = u^2 + c   
    
    @threads :greedy 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

julia_matrix_p3! (generic function with 1 method)

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]])

PyObject <matplotlib.image.AxesImage object at 0x7fe660306de0>

@benchmark julia_matrix_p3!(-0.8 + 0.256im, 10.0, $res, $ims, $data)

BenchmarkTools.Trial: 83 samples with 1 evaluation.
 Range (min … max):  50.895 ms … 83.096 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     58.658 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   60.496 ms ±  8.258 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

    █   ▂                       ▃                              
  ███████▁▄▄▁▁▅▁▅▅▄▅▅▄▄▄▄▁▄▄▁▅▇▁█▁▁▄▁▁▁▅▅▄▇▅▅▄▄▁▁▄▁▄▁▁▁▁▁▁▁▁▄ ▁
  50.9 ms         Histogram: frequency by time        80.5 ms <

 Memory estimate: 26.41 KiB, allocs estimate: 1461.

Zajímavé, zkusme změnit parametry.

@benchmark julia_matrix_p1!(-0.8 + 0.256im, 8.0, $res, $ims, $data, imax=5_000)

BenchmarkTools.Trial: 66 samples with 1 evaluation.
 Range (min … max):  61.300 ms … 117.568 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     77.108 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   76.921 ms ±  10.182 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

                           ▂▂█                                  
  ▄█▄▆▆█▄▄▁▄▁▁▁▁▁▁▄▄▄▄▄▁▄▁▆█████▆▆▄▆█▄▁▁▁▄▁▄▄▁▄▁▁▁▁▁▄▁▁▁▁▁▆▄▄▄ ▁
  61.3 ms         Histogram: frequency by time         96.3 ms <

 Memory estimate: 2.47 KiB, allocs estimate: 22.

@benchmark julia_matrix_p2!(-0.8 + 0.256im, 8.0, $res, $ims, $data, imax=5_000)

BenchmarkTools.Trial: 58 samples with 1 evaluation.
 Range (min … max):  72.189 ms … 108.980 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     85.155 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   86.703 ms ±   9.849 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁    ▁ █ ▄▁▄▄        ▁█     ▁▁  ▁     █▁  ▁           ▁       
  █▁▁▆▆█▁█▆████▁▆▆▆▁▁▆▁██▁▁▁▁▆██▁▁█▆▁▁▁▁██▆▁█▆▁▁▁▆▆▆▁▁▁▁█▆▁▁▁▆ ▁
  72.2 ms         Histogram: frequency by time          107 ms <

 Memory estimate: 2.47 KiB, allocs estimate: 22.

@benchmark julia_matrix_p3!(-0.8 + 0.256im, 8.0, $res, $ims, $data, imax=5_000)

BenchmarkTools.Trial: 72 samples with 1 evaluation.
 Range (min … max):  59.494 ms … 83.697 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     68.083 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   69.693 ms ±  5.469 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

               ▂▂     ▆█▆                           ▂          
  ▄▁▁▄▁▁▁▁▁▁▄▄███▆█▆█████▆▁▄▁▆▄▁▄▄▁▄▄▁▁▁▄▆▁▄▁▄▁▁▁▁▁▄█▆▄▁▁▁▁▁▄ ▁
  59.5 ms         Histogram: frequency by time        83.2 ms <

 Memory estimate: 26.41 KiB, allocs estimate: 1461.

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_p4!(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

julia_matrix_p4! (generic function with 1 method)

res = LinRange(-2, 2, 1920)
ims = LinRange(-1, 1, 1080)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p4!(-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]])

PyObject <matplotlib.image.AxesImage object at 0x7fe660156d80>

@benchmark julia_matrix_p4!(-0.8 + 0.256im, 8.0, $res, $ims, $data, imax=5_000)

BenchmarkTools.Trial: 85 samples with 1 evaluation.
 Range (min … max):  47.571 ms … 83.266 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     52.741 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   58.958 ms ± 10.872 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▃▄▃██▄                                    ▄▁              
  ▇▄▆██████▆▄▇▆▆▁▁▆▁▄▁▄▁▁▁▁▁▁▄▄▁▁▁▁▁▁▁▁▆▄▁▁▁▆▆▁██▄▄▄▆▄▁▁▁▁▁▁▄ ▁
  47.6 ms         Histogram: frequency by time        82.3 ms <

 Memory estimate: 2.80 KiB, allocs estimate: 29.

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_p4!(-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]])

PyObject <matplotlib.image.AxesImage object at 0x7fe6601ac680>

res = LinRange(-0.4, 0.4, 2000)
ims = LinRange(-0.6, -0.1, 1600)
data = zeros(Int64, length(res), length(ims))
julia_matrix_p4!(-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]])

PyObject <matplotlib.image.AxesImage object at 0x7fe65fd2bd70>

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ýraz EXPRESSION na procesu s id ID (může být rovno :any, pokud nám nezáleží na kterém), vrací objekt typu Future.
fetch(future), wait(future): získá výsledek vzdáleného výpočtu future, resp. počká na doběhnutí výpočtu.
@everywhere EXPRESSION: vyhodnotí výraz EXPRESSION 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átorem OP.

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()

workers()

4-element Vector{Int64}:
 2
 3
 4
 5

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 $0$ 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 $n=1\,000\,000$ krocích? Jak bude vypadat hustota pravděpodobnosti, že se po těchto $n$ krocích nachází v jistém bodě $k \in \mathbb{Z}$ ?

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 (generic function with 1 method)

random_walk(0)

random_walk(10)

-2

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)

10-element Vector{Int64}:
  14
  -6
  10
  42
  50
 -20
 -46
 -44
 -28
 -62

Kolik času nám zabere $10\,000$ procházek o stejném počtu kroků?

@time get_walks(10_000, 1_000_000);

 18.446609 seconds (3 allocations: 78.188 KiB)

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

Future(1, 1, 6, ReentrantLock(nothing, 0x00000000, 0x00, Base.GenericCondition{Base.Threads.SpinLock}(Base.IntrusiveLinkedList{Task}(nothing, nothing), Base.Threads.SpinLock(0)), (1, 0, 432325345510)), nothing)

fetch(t)

-4

t = @spawnat 2 random_walk(100) # v tento okamžik selže

Future(2, 1, 7, ReentrantLock(nothing, 0x00000000, 0x00, Base.GenericCondition{Base.Threads.SpinLock}(Base.IntrusiveLinkedList{Task}(nothing, nothing), Base.Threads.SpinLock(0)), (137176591061136, 137176591061168, 419807892823)), nothing)

fetch(t)

On worker 2:
UndefVarError: `#random_walk` not defined in `Main`
Suggestion: add an appropriate import or assignment. This global was declared but not assigned.
Stacktrace:
  [1] deserialize_datatype
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:1471
  [2] handle_deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:897
  [3] deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:844
  [4] handle_deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:904
  [5] deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:844 [inlined]
  [6] deserialize_global_from_main
    @ /usr/share/julia/stdlib/v1.11/Distributed/src/clusterserialize.jl:160
  [7] #5
    @ /usr/share/julia/stdlib/v1.11/Distributed/src/clusterserialize.jl:72 [inlined]
  [8] foreach
    @ ./abstractarray.jl:3187
  [9] deserialize
    @ /usr/share/julia/stdlib/v1.11/Distributed/src/clusterserialize.jl:72
 [10] handle_deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:990
 [11] deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:844
 [12] handle_deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:901
 [13] deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:844
 [14] handle_deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:904
 [15] deserialize
    @ /usr/share/julia/stdlib/v1.11/Serialization/src/Serialization.jl:844 [inlined]
 [16] deserialize_msg
    @ /usr/share/julia/stdlib/v1.11/Distributed/src/messages.jl:87
 [17] #invokelatest#2
    @ ./essentials.jl:1055 [inlined]
 [18] invokelatest
    @ ./essentials.jl:1052 [inlined]
 [19] message_handler_loop
    @ /usr/share/julia/stdlib/v1.11/Distributed/src/process_messages.jl:176
 [20] process_tcp_streams
    @ /usr/share/julia/stdlib/v1.11/Distributed/src/process_messages.jl:133
 [21] #103
    @ /usr/share/julia/stdlib/v1.11/Distributed/src/process_messages.jl:121

Stacktrace:
 [1] remotecall_fetch(f::Function, w::Distributed.Worker, args::Distributed.RRID; kwargs::@Kwargs{})
   @ Distributed /usr/share/julia/stdlib/v1.11/Distributed/src/remotecall.jl:465
 [2] remotecall_fetch(f::Function, w::Distributed.Worker, args::Distributed.RRID)
   @ Distributed /usr/share/julia/stdlib/v1.11/Distributed/src/remotecall.jl:454
 [3] remotecall_fetch
   @ /usr/share/julia/stdlib/v1.11/Distributed/src/remotecall.jl:492 [inlined]
 [4] call_on_owner
   @ /usr/share/julia/stdlib/v1.11/Distributed/src/remotecall.jl:565 [inlined]
 [5] fetch(r::Future)
   @ Distributed /usr/share/julia/stdlib/v1.11/Distributed/src/remotecall.jl:619
 [6] top-level scope
   @ In[27]:1

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)

Future(3, 1, 19, ReentrantLock(nothing, 0x00000000, 0x00, Base.GenericCondition{Base.Threads.SpinLock}(Base.IntrusiveLinkedList{Task}(nothing, nothing), Base.Threads.SpinLock(0)), (0, 0, 5)), nothing)

fetch(t)

-18

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)

10-element Vector{Int64}:
 2
 3
 4
 5
 3
 2
 3
 2
 3
 2

workers()

4-element Vector{Int64}:
 2
 3
 4
 5

myid()

function get_walks_p1(n, steps)
    output = fill(steps, n)
    output = pmap(random_walk, output)
    return output
end

get_walks_p1 (generic function with 1 method)

get_walks_p1(10, 1_000)

10-element Vector{Int64}:
  -2
 -26
  -4
   4
  14
  28
 -34
  -4
 -20
   2

Kolik času nám zabere $10\,000$ procházek o stejném počtu kroků?

@time get_walks_p1(10_000, 1_000_000);

  6.907404 seconds (846.99 k allocations: 41.633 MiB, 0.36% gc time, 2 lock conflicts)

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 (generic function with 1 method)

get_walks_p2(10, 1_000)

10-element SharedVector{Int64}:
  48
 -20
   8
  -6
   6
  40
 -36
 -34
  34
  26

Kolik času nám zabere $10\,000$ procházek o stejném počtu kroků?

@time get_walks_p2(10_000, 1_000_000);

  5.276671 seconds (1.11 k allocations: 58.148 KiB)

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í $f(n) = \begin{cases} n/2, & n \text{ je sudé}, \\ 3n+1, & n \text{ je liché}. \end{cases}$ Collatzova hypotéza pak zní: Pro každé přirozené $n$ existuje přirozené $k$ takové, že $f^{\circ k}(n) = 1$ .

nworkers()

Definujme metodu collatz, která pro dané $n$ provede Collatzovu iteraci a vrátí počet iterací nutný k dosažení $1$

@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 $p(x) = \sum_{k=0}^n a_k x^k$ 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.14069878366063912

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
=========================================================
  31╎31   @Base/math.jl:0; pow_body(x::Float64, n::Int64)
    ╎8394 …Julia/src/eventloop.jl:38; (::IJulia.var"#15#18")()
    ╎ 8394 …Julia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
    ╎  8394 @Base/essentials.jl:1052; invokelatest
    ╎   8394 @Base/essentials.jl:1055; #invokelatest#2
    ╎    8394 …c/execute_request.jl:74; execute_request(socket::ZMQ.Socket, msg…
    ╎     8394 …c/SoftGlobalScope.jl:65; softscope_include_string(m::Module, co…
    ╎    ╎ 8394 @Base/loading.jl:2643; include_string(mapexpr::typeof(REPL.soft…
   6╎    ╎  8394 @Base/boot.jl:430; eval
    ╎    ╎   8250 In[44]:4; poly_eval_v1(coeffs::Vector{Float64}, x::Float64)
  27╎    ╎    27   …ase/essentials.jl:916; getindex
  36╎    ╎    36   …ase/essentials.jl:917; getindex
 111╎    ╎    111  @Base/float.jl:493; *
   2╎    ╎    2    @Base/float.jl:491; +
    ╎    ╎    8074 @Base/math.jl:1199; ^
 178╎    ╎     178  @Base/math.jl:0; pow_body(x::Float64, n::Int64)
  36╎    ╎     36   @Base/math.jl:1202; pow_body(x::Float64, n::Int64)
    ╎    ╎     3    @Base/math.jl:1205; pow_body(x::Float64, n::Int64)
   3╎    ╎    ╎ 3    @Base/promotion.jl:639; ==
 334╎    ╎     633  @Base/math.jl:1213; pow_body(x::Float64, n::Int64)
    ╎    ╎    ╎ 299  @Base/operators.jl:379; >
 299╎    ╎    ╎  299  @Base/int.jl:83; <
    ╎    ╎     187  @Base/math.jl:1214; pow_body(x::Float64, n::Int64)
    ╎    ╎    ╎ 187  @Base/operators.jl:379; >
 187╎    ╎    ╎  187  @Base/int.jl:83; <
    ╎    ╎     851  @Base/math.jl:1215; pow_body(x::Float64, n::Int64)
  52╎    ╎    ╎ 52   @Base/float.jl:493; *
 798╎    ╎    ╎ 799  @Base/float.jl:496; muladd
    ╎    ╎     1456 @Base/math.jl:1216; pow_body(x::Float64, n::Int64)
    ╎    ╎    ╎ 1219 @Base/math.jl:56; two_mul
1219╎    ╎    ╎  1219 @Base/float.jl:493; *
    ╎    ╎    ╎ 237  @Base/math.jl:57; two_mul
    ╎    ╎    ╎  237  …ase/floatfuncs.jl:357; fma
 237╎    ╎    ╎   237  …se/floatfuncs.jl:352; fma_llvm
    ╎    ╎     1370 @Base/math.jl:1217; pow_body(x::Float64, n::Int64)
1370╎    ╎    ╎ 1370 @Base/float.jl:491; +
    ╎    ╎     636  @Base/math.jl:1219; pow_body(x::Float64, n::Int64)
    ╎    ╎    ╎ 636  @Base/operators.jl:596; *
 558╎    ╎    ╎  558  @Base/float.jl:493; *
    ╎    ╎    ╎  78   @Base/promotion.jl:430; *
  78╎    ╎    ╎   78   @Base/float.jl:493; *
    ╎    ╎     1069 @Base/math.jl:1220; pow_body(x::Float64, n::Int64)
    ╎    ╎    ╎ 965  @Base/math.jl:56; two_mul
 965╎    ╎    ╎  965  @Base/float.jl:493; *
    ╎    ╎    ╎ 104  @Base/math.jl:57; two_mul
    ╎    ╎    ╎  104  …ase/floatfuncs.jl:357; fma
 104╎    ╎    ╎   104  …se/floatfuncs.jl:352; fma_llvm
    ╎    ╎     96   @Base/math.jl:1221; pow_body(x::Float64, n::Int64)
  96╎    ╎    ╎ 96   @Base/float.jl:491; +
    ╎    ╎     1145 @Base/math.jl:1222; pow_body(x::Float64, n::Int64)
    ╎    ╎    ╎ 1145 @Base/int.jl:538; >>>
1145╎    ╎    ╎  1145 @Base/int.jl:530; >>>
    ╎    ╎     12   @Base/math.jl:1224; pow_body(x::Float64, n::Int64)
   1╎    ╎    ╎ 1    @Base/float.jl:493; *
  11╎    ╎    ╎ 11   @Base/float.jl:496; muladd
    ╎    ╎     402  @Base/math.jl:1225; pow_body(x::Float64, n::Int64)
 301╎    ╎    ╎ 301  …ase/essentials.jl:796; ifelse
    ╎    ╎    ╎ 101  @Base/float.jl:705; isfinite
 101╎    ╎    ╎  101  @Base/float.jl:492; -
    ╎    ╎   105  In[44]:5; poly_eval_v1(coeffs::Vector{Float64}, x::Float64)
    ╎    ╎    105  @Base/range.jl:908; iterate
 105╎    ╎     105  @Base/promotion.jl:639; ==
   1╎    ╎   1    @Base/math.jl:0; pow_body(x::Float64, n::Int64)
  32╎    ╎   32   @Base/math.jl:1202; pow_body(x::Float64, n::Int64)
   1╎1    …Julia/src/heartbeat.jl:22; heartbeat_thread(sock::Ptr{Nothing})
Total snapshots: 8426. 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
 =====  ======== ====                    ==== ========
  8250         0 In[44]                     4 poly_eval_v1(coeffs::Vector{Float…
   105         0 In[44]                     5 poly_eval_v1(coeffs::Vector{Float…
  8394         6 @Base/boot.jl            430 eval
  8394         0 @Base/essentials.jl     1055 #invokelatest#2
    27        27 @Base/essentials.jl      916 getindex
    36        36 @Base/essentials.jl      917 getindex
   301       301 @Base/essentials.jl      796 ifelse
  8394         0 @Base/essentials.jl     1052 invokelatest
  2984      2984 @Base/float.jl           493 *
  1468      1468 @Base/float.jl           491 +
   101       101 @Base/float.jl           492 -
   101         0 @Base/float.jl           705 isfinite
   810       809 @Base/float.jl           496 muladd
   341         0 @Base/floatfuncs.jl      357 fma
   341       341 @Base/floatfuncs.jl      352 fma_llvm
   486       486 @Base/int.jl              83 <
  1145      1145 @Base/int.jl             530 >>>
  1145         0 @Base/int.jl             538 >>>
  8394         0 @Base/loading.jl        2643 include_string(mapexpr::typeof(RE…
  8074         0 @Base/math.jl           1199 ^
   210       210 @Base/math.jl              ? pow_body(x::Float64, n::Int64)
    68        68 @Base/math.jl           1202 pow_body(x::Float64, n::Int64)
     3         0 @Base/math.jl           1205 pow_body(x::Float64, n::Int64)
   633       334 @Base/math.jl           1213 pow_body(x::Float64, n::Int64)
   187         0 @Base/math.jl           1214 pow_body(x::Float64, n::Int64)
   851         0 @Base/math.jl           1215 pow_body(x::Float64, n::Int64)
  1456         0 @Base/math.jl           1216 pow_body(x::Float64, n::Int64)
  1370         0 @Base/math.jl           1217 pow_body(x::Float64, n::Int64)
   636         0 @Base/math.jl           1219 pow_body(x::Float64, n::Int64)
  1069         0 @Base/math.jl           1220 pow_body(x::Float64, n::Int64)
    96         0 @Base/math.jl           1221 pow_body(x::Float64, n::Int64)
  1145         0 @Base/math.jl           1222 pow_body(x::Float64, n::Int64)
    12         0 @Base/math.jl           1224 pow_body(x::Float64, n::Int64)
   402         0 @Base/math.jl           1225 pow_body(x::Float64, n::Int64)
  2184         0 @Base/math.jl             56 two_mul
   341         0 @Base/math.jl             57 two_mul
   636         0 @Base/operators.jl       596 *
   486         0 @Base/operators.jl       379 >
    78         0 @Base/promotion.jl       430 *
   108       108 @Base/promotion.jl       639 ==
   105         0 @Base/range.jl           908 iterate
  8394         0 …ulia/src/eventloop.jl    38 (::IJulia.var"#15#18")()
  8394         0 …ulia/src/eventloop.jl     8 eventloop(socket::ZMQ.Socket)
  8394         0 …rc/execute_request.jl    74 execute_request(socket::ZMQ.Socke…
     1         1 …ulia/src/heartbeat.jl    22 heartbeat_thread(sock::Ptr{Nothin…
  8394         0 …rc/SoftGlobalScope.jl    65 softscope_include_string(m::Modul…
Total snapshots: 8426. 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

@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.14069878366063895

Profile.print()

Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
    ╎8438 …Julia/src/eventloop.jl:38; (::IJulia.var"#15#18")()
    ╎ 8438 …Julia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
    ╎  8438 @Base/essentials.jl:1052; invokelatest
    ╎   8438 @Base/essentials.jl:1055; #invokelatest#2
    ╎    8438 …c/execute_request.jl:74; execute_request(socket::ZMQ.Socket, msg…
    ╎     8438 …c/SoftGlobalScope.jl:65; softscope_include_string(m::Module, co…
    ╎    ╎ 8438 @Base/loading.jl:2643; include_string(mapexpr::typeof(REPL.soft…
    ╎    ╎  8438 @Base/boot.jl:430; eval
    ╎    ╎   4172 In[54]:6; poly_eval_v2(coeffs::Vector{Float64}, x::Float64)
  62╎    ╎    62   @Base/essentials.jl:916; getindex
  39╎    ╎    39   @Base/float.jl:493; *
4070╎    ╎    4071 @Base/float.jl:491; +
    ╎    ╎   320  In[54]:7; poly_eval_v2(coeffs::Vector{Float64}, x::Float64)
 320╎    ╎    320  @Base/float.jl:493; *
    ╎    ╎   3946 In[54]:8; poly_eval_v2(coeffs::Vector{Float64}, x::Float64)
    ╎    ╎    3946 @Base/range.jl:908; iterate
3946╎    ╎     3946 @Base/promotion.jl:639; ==
Total snapshots: 8438. 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 $a x^3 + b x^2 + c x + d = \big((ax + b)x + c\big)x + d.$

"""
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.14069878366063915

Profile.print()

Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
   ╎267 @IJulia/src/eventloop.jl:38; (::IJulia.var"#15#18")()
   ╎ 267 @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
   ╎  267 @Base/essentials.jl:1052; invokelatest
   ╎   267 @Base/essentials.jl:1055; #invokelatest#2
   ╎    267 …rc/execute_request.jl:74; execute_request(socket::ZMQ.Socket, msg:…
   ╎     267 …rc/SoftGlobalScope.jl:65; softscope_include_string(m::Module, cod…
   ╎    ╎ 267 @Base/loading.jl:2643; include_string(mapexpr::typeof(REPL.softsc…
   ╎    ╎  267 @Base/boot.jl:430; eval
   ╎    ╎   101 In[59]:7; poly_eval_v3(coeffs::Vector{Float64}, x::Float64)
  1╎    ╎    1   @Base/essentials.jl:916; getindex
100╎    ╎    100 @Base/float.jl:491; +
   ╎    ╎   166 In[59]:8; poly_eval_v3(coeffs::Vector{Float64}, x::Float64)
   ╎    ╎    166 @Base/range.jl:908; iterate
165╎    ╎     166 @Base/promotion.jl:639; ==
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.14069878366063915

Profile.print()

Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
   ╎253 @IJulia/src/eventloop.jl:38; (::IJulia.var"#15#18")()
   ╎ 253 @IJulia/src/eventloop.jl:8; eventloop(socket::ZMQ.Socket)
   ╎  253 @Base/essentials.jl:1052; invokelatest
   ╎   253 @Base/essentials.jl:1055; #invokelatest#2
   ╎    253 …rc/execute_request.jl:74; execute_request(socket::ZMQ.Socket, msg:…
   ╎     253 …rc/SoftGlobalScope.jl:65; softscope_include_string(m::Module, cod…
   ╎    ╎ 253 @Base/loading.jl:2643; include_string(mapexpr::typeof(REPL.softsc…
   ╎    ╎  253 @Base/boot.jl:430; eval
   ╎    ╎   224 In[66]:7; poly_eval_v4(reversed_coeffs::Vector{Float64}, x::Flo…
106╎    ╎    106 @Base/float.jl:493; *
118╎    ╎    118 @Base/float.jl:491; +
 28╎    ╎   29  In[66]:8; poly_eval_v4(reversed_coeffs::Vector{Float64}, x::Flo…
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: 31 samples with 1 evaluation.
 Range (min … max):  154.850 ms … 168.707 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     160.125 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   160.709 ms ±   3.426 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

                ▃ ▃    ▃ █ ▃                                     
  ▇▁▁▇▁▇▁▇▁▁▇▁▁▁█▁█▁▇▇▇█▁█▇█▁▇▁▇▁▁▁▇▇▁▇▁▁▁▁▇▁▁▇▇▁▁▁▇▁▁▁▁▁▁▁▁▇▁▇ ▁
  155 ms           Histogram: frequency by time          169 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark poly_eval_v2($coeffs, -0.5)

BenchmarkTools.Trial: 3104 samples with 1 evaluation.
 Range (min … max):  1.411 ms …   2.842 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.531 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.600 ms ± 195.952 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅█▃ ▁▃                                                       
  ███████▇▆▆▆▅▅▅▄▄▄▃▃▃▃▃▄▃▃▃▂▂▂▂▂▂▂▂▂▂▁▂▁▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  1.41 ms         Histogram: frequency by time         2.3 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark poly_eval_v3($coeffs, -0.5)

BenchmarkTools.Trial: 1768 samples with 1 evaluation.
 Range (min … max):  2.408 ms …   5.177 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.623 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.812 ms ± 467.900 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▂▂▁                                                         
  █████▇▇▄▅▅▄▃▃▃▃▂▂▃▃▂▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▁▁▂▂▁▂▂▁▁▂▁▁▁▁▁▁▁▁▁▁ ▂
  2.41 ms         Histogram: frequency by time        4.36 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark poly_eval_v4($coeffs, -0.5)

BenchmarkTools.Trial: 1876 samples with 1 evaluation.
 Range (min … max):  2.387 ms …   4.893 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.530 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.654 ms ± 333.582 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █ ▆ ▁                                                        
  ███████▆▄▄▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  2.39 ms         Histogram: frequency by time        3.96 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: 3648 samples with 1 evaluation.
 Range (min … max):  1.206 ms …   2.205 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.296 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.361 ms ± 168.389 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃█                                                           
  ███▄▆▆▆▅▄▄▄▃▃▃▃▃▃▃▂▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  1.21 ms         Histogram: frequency by time        1.92 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: 1913 samples with 1 evaluation.
 Range (min … max):  2.390 ms …   3.537 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.537 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.602 ms ± 202.174 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▆█▁ ▃▂▁                                                     
  ▄███████▇██▇▇▆▇▆▅▅▅▄▄▅▄▄▃▃▄▃▃▃▃▃▃▂▂▃▃▃▃▃▃▃▃▂▃▂▂▃▂▃▂▂▂▂▂▂▂▁▂ ▃
  2.39 ms         Histogram: frequency by time        3.28 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.

10: Paralelní výpočty a profilování

1. Úvod

2. Paralelní výpočty

2.1: Více vláknové výpočty

Threads.@threads makro

Threads.@spawn makro, fetch a wait

Cvičení, resp. poznámka k 1. domácímu úkolu (z roku 2021)

Cvičení: Paralelizace výpočtu Julia fraktálu

2.2 Distribuované výpočty

Demonstrační ukázka: náhodná procházka

Cvičení: Collatz conjecture

3. Profilování

3.1 Ukázka

3.2 Vylepšení I

3.3 Vylepšení II

3.4 Poznámka

Řešení některých úkolů

Reference

`Threads.@threads` makro

`Threads.@spawn` makro, `fetch` a `wait`