Add lecture 10

Author: VaclavMacha

docs_vitepress/make.jl

@@ -75,11 +75,11 @@ pages = [
         #     "Lecture" => "lecture.md",
         #     "Lab" => "lab.md",
         # ]),
-        # "10: Parallel programming" => add_prefix("lecture_10", [
-        #     "Lecture" => "lecture.md",
-        #     "Lab" => "lab.md",
-        #     "Homework" => "hw.md",
-        # ]),
+        "10: Parallel programming" => add_prefix("lecture_10", [
+            "Lecture" => "lecture.md",
+            "Lab" => "lab.md",
+            "Homework" => "hw.md",
+        ]),
         "11: GPU programming" => add_prefix("lecture_11", [
             "Lecture" => "lecture.md",
             "Lab" => "lab.md",

docs_vitepress/src/lectures/lecture_10/42-years-processor-trend.png

@@ -0,0 +1,75 @@
+@everywhere begin 
+	"""
+		sample_all_installed_pkgs(path::AbstractString)
+
+	Returns root folders of all installed packages in the system. Package version is sampled.
+	"""
+	function sample_all_installed_pkgs(path::AbstractString)
+		pkgs = readdir(path)
+		# [rand(readdir(joinpath(path, p), join=true)) for p in pkgs] # sampling version
+		[readdir(joinpath(path, p), join=true)[1] for p in pkgs if isdir(joinpath(path, p))]    # deterministic version
+	end
+
+	"""
+		filter_jl(path)
+
+	Recursively walks the directory structure to obtain all `.jl` files.
+	"""
+	filter_jl(path) = reduce(vcat, joinpath.(rootpath, filter(endswith(".jl"), files)) for (rootpath, dirs, files) in walkdir(path))
+
+	"""
+		tokenize(jl_path)
+
+	Parses a ".jl" file located at `jl_path` and extracts all symbols and expression heads from the extracted AST.
+	"""
+	function tokenize(jl_path)
+		_extract_symbols(x) = Symbol[]
+		_extract_symbols(x::Symbol) = [x]
+		function _extract_symbols(x::Expr) 
+			if length(x.args) > 0
+				Symbol.(vcat(x.head, reduce(vcat, _extract_symbols(arg) for arg in x.args)))
+			else
+				Symbol[]
+			end
+		end
+		
+		scode = "begin\n" * read(jl_path, String) * "end\n"
+		try 
+			code = Meta.parse(scode)
+			_extract_symbols(code)
+		catch e
+			if ~isa(e, Meta.ParseError)
+				rethrow(e)		
+			end
+			Symbol[]
+		end
+	end
+
+
+	function histtokens!(h, filename::AbstractString)
+		for t in tokenize(filename)
+			h[t] = get(h, t, 0) + 1
+		end
+		h
+	end
+
+	function dohistogram(chnl)
+		h = Dict{Symbol, Int}()
+		while isready(chnl)
+			f = take!(chnl)
+			histtokens!(h, f)
+		end
+		return(h)
+	end
+end
+
+chnl = RemoteChannel() do 
+	Channel(typemax(Int)) do ch
+		for package in sample_all_installed_pkgs("/Users/tomas.pevny/.julia/packages")
+		    foreach(c -> put!(ch, c), filter_jl(package))
+		end
+	end
+end
+
+mapreduce(fetch, mergewith(+), [@spawnat i dohistogram(chnl) for i in workers()])
+

docs_vitepress/src/lectures/lecture_10/exercise.jl

@@ -0,0 +1,43 @@
+# [Homework 9: Accelerating 1D convolution with threads](@id hw09)
+
+## How to submit
+
+Put all the code of inside `hw.jl`. Zip only this file (not its parent folder) and upload it to BRUTE. You should not not import anything but `Base.Threads` or just `Threads`. 
+
+::: danger Homework (2 points)
+
+Implement *multithreaded* discrete 1D convolution operator[^1] without padding (output will be shorter). The required function signature: `thread_conv1d(x, w)`, where `x` is the signal array and `w` the kernel. For testing correctness of the implementation you can use the following example of a step function and it's derivative realized by kernel `[-1, 1]`:
+
+```julia
+using Test
+@test all(thread_conv1d(vcat([0.0, 0.0, 1.0, 1.0, 0.0, 0.0]), [-1.0, 1.0]) .≈ [0.0, -1.0, 0.0, 1.0, 0.0])
+```
+
+[^1]: Discrete convolution with finite support [https://en.wikipedia.org/wiki/Convolution#Discrete\_convolution](https://en.wikipedia.org/wiki/Convolution#Discrete_convolution)
+
+Your parallel implementation will be tested both in sequential and two threaded mode with the following inputs
+
+```julia
+using Random
+Random.seed!(42)
+x = rand(10_000_000)
+w = [1.0, 2.0, 4.0, 2.0, 1.0]
+@btime thread_conv1d($x, $w);
+```
+
+On your local machine you should be able to achieve `0.6x` reduction in execution time with two threads, however the automatic eval system is a noisy environment and therefore we require only `0.8x` reduction therein. This being said, please reach out to us, if you encounter any issues.
+
+**HINTS**:
+
+- start with single threaded implementation
+- don't forget to reverse the kernel
+- `@threads` macro should be all you need
+- for testing purposes create a simple script, that you can run with `julia -t 1` and `julia -t 2`
+
+:::
+
+::: details Show solution
+
+Nothing to see here
+
+:::

docs_vitepress/src/lectures/lecture_10/hw.md

@@ -0,0 +1,75 @@
+using Pkg
+Pkg.activate(@__DIR__)
+using GLMakie
+using BenchmarkTools
+using Distributed
+using SharedArrays
+
+function juliaset_pixel(z₀, c)
+    z = z₀
+    for i in 1:255
+        abs2(z)> 4.0 && return (i - 1)%UInt8
+        z = z*z + c
+    end
+    return UInt8(255)
+end
+
+function juliaset_column!(img, c, n, colj, j)
+    x = -2.0 + (j-1)*4.0/(n-1)
+    for i in 1:n
+        y = -2.0 + (i-1)*4.0/(n-1)
+        @inbounds img[i, colj] = juliaset_pixel(x+im*y, c)
+    end
+    nothing
+end
+
+function juliaset_columns(c, n, columns)
+    img = Array{UInt8,2}(undef, n, length(columns))
+    for (colj, j) in enumerate(columns)
+        juliaset_column!(img, c, n, colj, j)
+    end
+    img
+end
+
+function juliaset_distributed(x, y, partitions = nworkers(), n = 1000)
+    c = x + y*im
+    columns = Iterators.partition(1:n, div(n, partitions))
+    slices = pmap(cols -> juliaset_columns(c, n, cols), columns)
+    reduce(hcat, slices)
+end
+
+# @btime juliaset_distributed(-0.79, 0.15)
+
+# frac = juliaset_distributed(-0.79, 0.15)
+# plot(heatmap(1:size(frac,1), 1:size(frac,2), frac, color=:Spectral))
+
+
+####
+#   Let's work out the shared array approach
+####
+function juliaset_column!(img, c, n, j)
+    x = -2.0 + (j-1)*4.0/(n-1)
+    for i in 1:n
+        y = -2.0 + (i-1)*4.0/(n-1)
+        @inbounds img[i, j] = juliaset_pixel(x+im*y, c)
+    end
+    nothing
+end
+
+function juliaset_range!(img, c, n, columns)
+    for j in columns
+        juliaset_column!(img, c, n, j)
+    end
+    nothing
+end
+
+function juliaset_shared(x, y, partitions = nworkers(), n = 1000)
+    c = x + y*im
+    columns = Iterators.partition(1:n, div(n, partitions))
+    img = SharedArray{UInt8,2}((n, n))
+    slices = pmap(cols -> juliaset_range!(img, c, n, cols), columns)
+    img
+end
+
+# juliaset_shared(-0.79, 0.15)
+# juliaset_shared(-0.79, 0.15, 16)

docs_vitepress/src/lectures/lecture_10/juliaset_p.jl

@@ -0,0 +1,77 @@
+using Plots, BenchmarkTools
+Threads.nthreads()
+function juliaset_pixel(z₀, c)
+    z = z₀
+    for i in 1:255
+        abs2(z)> 4.0 && return (i - 1)%UInt8
+        z = z*z + c
+    end
+    return UInt8(255)
+end
+
+function juliaset_column!(img, c, n, j)
+    x = -2.0 + (j-1)*4.0/(n-1)
+    for i in 1:n
+        y = -2.0 + (i-1)*4.0/(n-1)
+        @inbounds img[i,j] = juliaset_pixel(x+im*y, c)
+    end
+    nothing
+end
+
+function juliaset_single!(img, c, n)
+    for j in 1:n
+        juliaset_column!(img, c, n, j)
+    end
+    nothing
+end
+
+function juliaset(x, y, n=1000, method = juliaset_single!, extra...)
+    c = x + y*im
+    img = Array{UInt8,2}(undef,n,n)
+    method(img, c, n, extra...)
+    return img
+end
+
+# frac = juliaset(-0.79, 0.15)
+# plot(heatmap(1:size(frac,1),1:size(frac,2), frac, color=:Spectral))
+
+
+@btime juliaset(-0.79, 0.15);
+
+function juliaset_static!(img, c, n)
+    Threads.@threads for j in 1:n
+        juliaset_column!(img, c, n, j)
+    end
+    nothing
+end
+
+@btime juliaset(-0.79, 0.15, 1000, juliaset_static!);
+
+
+using Folds
+function juliaset_folds!(img, c, n)
+    Folds.foreach(j -> juliaset_column!(img, c, n, j), 1:n)
+    nothing
+end
+julia> @btime juliaset(-0.79, 0.15, 1000, juliaset_folds!);
+  16.267 ms (25 allocations: 978.20 KiB)
+
+function juliaset_folds!(img, c, n, nt)
+    parts = collect(Iterators.partition(1:n, cld(n, nt)))
+    Folds.foreach(parts) do ii
+        foreach(j ->juliaset_column!(img, c, n, j), ii)
+    end
+    nothing
+end
+julia> @btime juliaset(-0.79, 0.15, 1000, (args...) -> juliaset_folds!(args..., 16));
+  16.716 ms (25 allocations: 978.61 KiB)
+
+
+using FLoops, FoldsThreads
+function juliaset_folds!(img, c, n)
+    @floop ThreadedEx(basesize = 2) for j in 1:n
+        juliaset_column!(img, c, n, j)
+    end
+    nothing
+end
+@btime juliaset(-0.79, 0.15, 1000, juliaset_folds!);