diff --git a/notebooks/LEQ.ipynb b/notebooks/LEQ.ipynb new file mode 100644 index 0000000..a7e4fca --- /dev/null +++ b/notebooks/LEQ.ipynb @@ -0,0 +1,163 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "599913a3", + "metadata": {}, + "source": [ + "# Solving Linear Equations\n", + "\n", + "## Serial Algorithm\n", + "First, we construct a linear equations system $Ax=b$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe3b5b02", + "metadata": {}, + "outputs": [], + "source": [ + "n = 5\n", + "A = rand(-10.0:10.0, (n, n))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57c912fb", + "metadata": {}, + "outputs": [], + "source": [ + "x = rand(-10.0:10.0, n)\n", + "b = A * x" + ] + }, + { + "cell_type": "markdown", + "id": "4baa0681", + "metadata": {}, + "source": [ + "The code in the following cell converts the problem $Ax=b$ to the upper triangular equation system $Ux=y$. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11d255e0", + "metadata": {}, + "outputs": [], + "source": [ + "function convert_to_upper_triangular(A,b)\n", + " # Upper Triangularization: convert Ax=b to Ux=y\n", + " for k in 1:n\n", + " for j in k+1:n\n", + " # Divide by pivot\n", + " A[k,j] = A[k,j] / A[k,k]\n", + " end\n", + " b[k] = b[k] / A[k,k]\n", + " A[k,k] = 1\n", + " # Substract lower rows\n", + " for i in k+1:n \n", + " for j in k+1:n\n", + " A[i,j]=A[i,j] - A[i,k] * A[k,j]\n", + " end\n", + " b[i] = b[i] - A[i,k] * b[k]\n", + " A[i,k] = 0\n", + " end\n", + " end\n", + " return A, b #U,y\n", + "end\n" + ] + }, + { + "cell_type": "markdown", + "id": "02c47593", + "metadata": {}, + "source": [ + "The function in the following cell solves the upper triangular equation system using backwards substitution. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3c19497c", + "metadata": {}, + "outputs": [], + "source": [ + "function solve_upper_triangular(U,y)\n", + " n = size(U,1)\n", + " for step in reverse(1:n)\n", + " if U[step,step] == 0\n", + " if y[step] != 0\n", + " return \"No solution\"\n", + " else\n", + " return \"Infinity solutions\"\n", + " end\n", + " else\n", + " # Backwards substitution\n", + " y[step] = y[step] / U[step,step]\n", + " end\n", + " for row in reverse(1:step-1)\n", + " y[row] -= U[row,step] * y[step]\n", + " end\n", + " end\n", + " return y \n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2c24d85e", + "metadata": {}, + "outputs": [], + "source": [ + "U,y = convert_to_upper_triangular(A,b)\n", + "sol = solve_upper_triangular(U,y)" + ] + }, + { + "cell_type": "markdown", + "id": "1c356b5a", + "metadata": {}, + "source": [ + "We can test if the obtained solution is correct using `@test`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c71828d", + "metadata": {}, + "outputs": [], + "source": [ + "using Test\n", + "@test sol ≈ x" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe3c5374", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.9.1", + "language": "julia", + "name": "julia-1.9" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.9.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/asp.ipynb b/notebooks/asp.ipynb index cedc122..843c863 100644 --- a/notebooks/asp.ipynb +++ b/notebooks/asp.ipynb @@ -31,7 +31,7 @@ "id": "ade31d26", "metadata": {}, "source": [ - "## The all pairs of shortest paths (ASP) problem\n", + "## The All Pairs of Shortest Paths (ASP) problem\n", "\n", "Let us start by presenting the all pairs of shortest paths (ASP) problem and its solution with the [Floyd–Warshall algorithm](https://en.wikipedia.org/wiki/Floyd%E2%80%93Warshall_algorithm).\n", "\n", @@ -70,10 +70,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "4fe447c5", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "floyd! (generic function with 1 method)" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "function floyd!(C)\n", " n = size(C,1)\n", @@ -99,10 +110,25 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "860e537c", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "4×4 Matrix{Int64}:\n", + " 0 9 6 1\n", + " 2 0 8 3\n", + " 5 3 0 6\n", + " 10 8 5 0" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "inf = 1000\n", "C = [\n", @@ -296,7 +322,7 @@ "id": "db9a3294", "metadata": {}, "source": [ - "TODO explain this part" + "At each iteration, the processor needs the input values `C[i,j]`, `C[i,k]` and `C[k,j]`. Since we split the data row-wise, the process already has values `C[i,j]` and `C[i,k]`. However, `C[k,j]` may be stored in a different process. " ] }, { @@ -314,6 +340,14 @@ "" ] }, + { + "cell_type": "markdown", + "id": "a01872ef", + "metadata": {}, + "source": [ + "As we iterate over columns $j$, the process needs input from all values of row $k$. Therefore, at the start of iteration $k$, the whole row $k$ needs to be communicated." + ] + }, { "attachments": { "g12957.png": { @@ -634,7 +668,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 1.9.0", + "display_name": "Julia 1.9.1", "language": "julia", "name": "julia-1.9" }, @@ -642,7 +676,7 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.9.0" + "version": "1.9.1" } }, "nbformat": 4, diff --git a/notebooks/solutions.ipynb b/notebooks/solutions.ipynb index 64b99de..779d2d7 100644 --- a/notebooks/solutions.ipynb +++ b/notebooks/solutions.ipynb @@ -291,6 +291,194 @@ "rmprocs(workers())" ] }, + { + "cell_type": "markdown", + "id": "19641daf", + "metadata": {}, + "source": [ + "## TSP Exercise: Measure search overhead" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f00557a0", + "metadata": {}, + "outputs": [], + "source": [ + "## TSP serial \n", + "function tsp_serial(connections,city)\n", + " num_cities = length(connections)\n", + " path=zeros(Int,num_cities)\n", + " hops = 1\n", + " path[hops] = city\n", + " min_path = zeros(Int, num_cities)\n", + " current_distance = 0\n", + " min_distance = typemax(Int)\n", + " # Collect search time \n", + " search_time = @elapsed min_path, min_distance = tsp_serial_impl(connections,hops,path,current_distance, min_path, min_distance)\n", + " (;path=min_path,distance=min_distance, search_time)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30784da2", + "metadata": {}, + "outputs": [], + "source": [ + "## TSP distributed\n", + "@everywhere function tsp_dist_impl(wait_time, connections,hops,path,current_distance,min_dist_chnl, max_hops,jobs_chnl,ftr_result)\n", + " num_cities = length(connections)\n", + " if hops == num_cities\n", + " min_distance = fetch(min_dist_chnl)\n", + " if current_distance < min_distance\n", + " take!(min_dist_chnl)\n", + " # Collect wait time to substract from overall search time \n", + " if ftr_result !== nothing\n", + " wait_time += @elapsed @spawnat 1 begin\n", + " result = fetch(ftr_result)\n", + " result.path .= path\n", + " result.min_distance_ref[] = current_distance\n", + " end |> wait\n", + " end\n", + " put!(min_dist_chnl, current_distance)\n", + " end\n", + " elseif hops <= max_hops\n", + " current_city = path[hops]\n", + " next_hops = hops + 1\n", + " for (next_city,distance_increment) in connections[current_city]\n", + " if !visited(next_city,hops,path)\n", + " path[next_hops] = next_city\n", + " next_distance = current_distance + distance_increment\n", + " # Collect wait time because fetch may block\n", + " wait_time += @elapsed min_distance = fetch(min_dist_chnl)\n", + " if next_distance < min_distance\n", + " tsp_dist_impl(wait_time, connections,next_hops,path,next_distance,min_dist_chnl,max_hops,jobs_chnl,ftr_result)\n", + " end\n", + " end\n", + " end \n", + " else\n", + " # Collect communication time and add to wait time\n", + " wait_time += @elapsed if jobs_chnl !== nothing \n", + " path_copy = copy(path) \n", + " put!(jobs_chnl,(;hops,path=path_copy,current_distance))\n", + " end\n", + " end\n", + " # Return wait time\n", + " wait_time\n", + "end\n", + "\n", + "function tsp_dist(connections,city)\n", + " max_hops = 2\n", + " num_cities = length(connections)\n", + " path=zeros(Int,num_cities)\n", + " result_path=zeros(Int, num_cities)\n", + " wait_time = 0\n", + " search_time = 0\n", + " hops = 1\n", + " path[hops] = city\n", + " current_distance = 0\n", + " min_distance = typemax(Int)\n", + " jobs_chnl = RemoteChannel(()->Channel{Any}(10))\n", + " min_dist_chnl = RemoteChannel(()->Channel{Int}(1))\n", + " put!(min_dist_chnl, min_distance)\n", + " ftr_result = @spawnat 1 (;path=result_path,min_distance_ref=Ref(min_distance))\n", + " @async begin\n", + " # Collect search time from master process\n", + " search_time += @elapsed wait_time += tsp_dist_impl(wait_time,connections,hops,path,current_distance,min_dist_chnl,max_hops,jobs_chnl,nothing)\n", + " for w in workers()\n", + " put!(jobs_chnl,nothing)\n", + " end\n", + " end\n", + " @sync for w in workers()\n", + " @spawnat w begin\n", + " path = zeros(Int, num_cities)\n", + " max_hops = typemax(Int)\n", + " while true\n", + " job = take!(jobs_chnl)\n", + " if job == nothing\n", + " break\n", + " end\n", + " hops = job.hops\n", + " path = job.path \n", + " current_distance = job.current_distance\n", + " min_distance = fetch(min_dist_chnl)\n", + " if current_distance < min_distance\n", + " # Collect search time from worker processes \n", + " search_time += @elapsed wait_time += tsp_dist_impl(wait_time,connections,hops,path,current_distance,min_dist_chnl,max_hops,nothing,ftr_result)\n", + " end\n", + " end\n", + " end\n", + " end \n", + " result = fetch(ftr_result)\n", + " (;path = result.path, distance = result.min_distance_ref[], search_time, wait_time)\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "694de934", + "metadata": {}, + "outputs": [], + "source": [ + "using Distributed\n", + "using RandomMatrix\n", + "using Plots\n", + "\n", + "function generate_rand_connections(city_range, distance_range)\n", + " # generate random connections matrix \n", + " n_cities = rand(city_range)\n", + " matrix = randTriangular(distance_range, n_cities; Diag=false)\n", + "\n", + " connections = Array{Array{Tuple{Int64,Int64},1},1}(undef, n_cities)\n", + " for i in 1:n_cities\n", + " connections[i] = Array{Tuple{Int64,Int64},1}(undef, n_cities)\n", + " end\n", + " for i in 1:n_cities\n", + " for j in i:n_cities\n", + " distance = matrix[i,j]\n", + " connections[i][j] = (j,distance)\n", + " connections[j][i] = (i,distance)\n", + " end\n", + " end\n", + " return connections\n", + "end\n", + "\n", + "# Run once so compile times are not measured\n", + "distance_range = 1:100\n", + "connections = generate_rand_connections(4:4, distance_range)\n", + "tsp_dist(connections,1)\n", + "tsp_serial(connections,1)\n", + "\n", + "# Measure runtimes of serial and parallel algorithm\n", + "n_it = 5\n", + "city_ranges = [4:4, 6:6, 8:8, 10:10]\n", + "search_overhead = zeros(Float64, length(city_ranges), n_it )\n", + "for (i, n) in enumerate(city_ranges)\n", + " for k in 1:n_it\n", + " connections = generate_rand_connections(n, distance_range)\n", + " @show n, k\n", + " path_dist, distance_dist, search_time_dist, wait_time_dist = tsp_dist(connections,1)\n", + " path_serial, distance_serial, search_time_serial = tsp_serial(connections,1)\n", + " # Compute search overhead as difference between distributed program and serial program\n", + " # (without time spent communicating or waiting)\n", + " search_overhead[i, k] = search_time_dist - wait_time_dist - search_time_serial\n", + " end\n", + "end\n", + "\n", + "min_search_oh = minimum(search_overhead, dims=2)\n", + "city_sizes = [4,6,8,10]\n", + "plot(city_sizes, min_search_oh, yaxis=:log, seriestype=:scatter,legend=false)\n", + "plot!(city_sizes, min_search_oh, yaxis=:log, legend=false)\n", + "\n", + "xlabel!(\"Number of cities\")\n", + "ylabel!(\"Search overhead (s)\")\n", + "title!(\"Minimum search overhead for different problem sizes\")" + ] + }, { "cell_type": "markdown", "id": "47d88e7a", @@ -314,7 +502,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 1.9.0", + "display_name": "Julia 1.9.1", "language": "julia", "name": "julia-1.9" }, @@ -322,7 +510,7 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.9.0" + "version": "1.9.1" } }, "nbformat": 4, diff --git a/notebooks/tsp-parallel.jl b/notebooks/tsp-parallel.jl new file mode 100644 index 0000000..2d02526 --- /dev/null +++ b/notebooks/tsp-parallel.jl @@ -0,0 +1,122 @@ +using Distributed + +if procs() == workers() + addprocs(4) +end + +@everywhere function visited(city,hops,path) + for i = 1:hops + if path[i] == city + return true + end + end + return false +end + +# solution [1, 4, 5, 2, 3, 6], distance = 222 +connections = [ + [(1,0),(4,39),(5,76), (6,78),(3,94),(2,97)], + [(2,0),(5,25),(4,58),(3,62),(1,97),(6,109)], + [(3,0),(6,58),(2,62),(4,68),(5,70),(1,94)], + [(4,0),(5,38),(1,39),(2,58),(3,68),(6,78)], + [(5,0),(2,25),(4,38),(3,70),(1,76),(6,104)], + [(6,0),(3,58),(1,78),(4,78),(5,104),(2,109)] +] + +# Shortest route with start 1: 1-3-2-4 (distance: 7) +con2 = [ + [(1,0), (2,2), (3,3), (4,4)], + [(2,0), (4,1), (1,2), (3,3)], + [(3,0), (1,3), (2,3), (4,10)], + [(4,0), (2,1), (1,4), (3,10)] +] + + +## TSP distributed +@everywhere function tsp_dist_impl(connections,hops,path,current_distance,min_dist_chnl, max_hops,jobs_chnl,ftr_result) + num_cities = length(connections) + if hops == num_cities + min_distance = fetch(min_dist_chnl) + if current_distance < min_distance + take!(min_dist_chnl) + # Wait until results are written to future + if ftr_result !== nothing + @spawnat 1 begin + result = fetch(ftr_result) + result.path .= path + result.min_distance_ref[] = current_distance + end |> wait + end + # Unblock waiting processes + put!(min_dist_chnl, current_distance) + end + elseif hops <= max_hops + current_city = path[hops] + next_hops = hops + 1 + for (next_city,distance_increment) in connections[current_city] + if !visited(next_city,hops,path) + path[next_hops] = next_city + next_distance = current_distance + distance_increment + min_distance = fetch(min_dist_chnl) + if next_distance < min_distance + tsp_dist_impl(connections,next_hops,path,next_distance,min_dist_chnl,max_hops,jobs_chnl,ftr_result) + end + end + end + else + if jobs_chnl !== nothing + # Allocate new memory so paths are not overwritten in queue + path_copy = copy(path) + put!(jobs_chnl,(;hops,path=path_copy,current_distance)) + end + end +end + +function tsp_dist(connections,city) + max_hops = 2 + num_cities = length(connections) + path=zeros(Int,num_cities) + result_path=zeros(Int, num_cities) + hops = 1 + path[hops] = city + current_distance = 0 + min_distance = typemax(Int) + jobs_chnl = RemoteChannel(()->Channel{Any}(10)) + # Initialize min distance channel with Intmax + min_dist_chnl = RemoteChannel(()->Channel{Int}(1)) + put!(min_dist_chnl, min_distance) + # Future to store overall result + ftr_result = @spawnat 1 (;path=result_path,min_distance_ref=Ref(min_distance)) + @async begin + tsp_dist_impl(connections,hops,path,current_distance,min_dist_chnl,max_hops,jobs_chnl,nothing) + for w in workers() + put!(jobs_chnl,nothing) + end + end + @sync for w in workers() + @spawnat w begin + path = zeros(Int, num_cities) + max_hops = typemax(Int) + jobs_channel = nothing + while true + job = take!(jobs_chnl) + if job == nothing + break + end + hops = job.hops + path = job.path + current_distance = job.current_distance + # Prune this job if current distance exeeds search threshold + min_distance = fetch(min_dist_chnl) + if current_distance < min_distance + tsp_dist_impl(connections,hops,path,current_distance,min_dist_chnl,max_hops,jobs_channel,ftr_result) + + end + end + end + end + result = fetch(ftr_result) + (;path = result.path, distance = result.min_distance_ref[]) +end +city = 1 +tsp_dist(connections,city) \ No newline at end of file diff --git a/notebooks/tsp-serial.jl b/notebooks/tsp-serial.jl new file mode 100644 index 0000000..196d64a --- /dev/null +++ b/notebooks/tsp-serial.jl @@ -0,0 +1,66 @@ +using Distributed + +connections = [ + [(1,0),(4,39),(5,76), (6,78),(3,94),(2,97)], + [(2,0),(5,25),(4,58),(3,62),(1,97),(6,109)], + [(3,0),(6,58),(2,62),(4,68),(5,70),(1,94)], + [(4,0),(5,38),(1,39),(2,58),(3,68),(6,78)], + [(5,0),(2,25),(4,38),(3,70),(1,76),(6,104)], + [(6,0),(3,58),(1,78),(4,78),(5,104),(2,109)] +] + +# Shortest route with start 1: 1-3-2-4 (distance: 7) +con2 = [ + [(1,0), (2,2), (3,3), (4,4)], + [(2,0), (4,1), (1,2), (3,3)], + [(3,0), (1,3), (2,3), (4,10)], + [(4,0), (2,1), (1,4), (3,10)] +] + +@everywhere function visited(city,hops,path) + for i = 1:hops + if path[i] == city + return true + end + end + return false +end + +## TSP serial +function tsp_serial_impl(connections,hops,path,current_distance, min_path, min_distance) + num_cities = length(connections) + if hops == num_cities + if current_distance < min_distance + min_path .= path + return min_path, current_distance + end + else + current_city = path[hops] + next_hops = hops + 1 + for (next_city,distance_increment) in connections[current_city] + if !visited(next_city,hops,path) + path[next_hops] = next_city + next_distance = current_distance + distance_increment + if next_distance < min_distance + min_path, min_distance = tsp_serial_impl(connections,next_hops,path,next_distance,min_path,min_distance) + end + end + end + end + return min_path, min_distance +end + +function tsp_serial(connections,city) + num_cities = length(connections) + path=zeros(Int,num_cities) + hops = 1 + path[hops] = city + min_path = zeros(Int, num_cities) + current_distance = 0 + min_distance = typemax(Int) + min_path, min_distance = tsp_serial_impl(connections,hops,path,current_distance, min_path, min_distance) + (;path=min_path,distance=min_distance) +end + +city = 1 +tsp_serial(connections,city)