{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e3666bcd",
   "metadata": {},
   "source": [
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/3/39/VU_logo.png/800px-VU_logo.png?20161029201021\" width=\"350\">\n",
    "\n",
    "### Programming large-scale parallel systems\n",
    "\n",
    "# Traveling sales person"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d49d8e7",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "In this notebook, we will learn\n",
    "\n",
    "- How to parallelize the solution of the traveling sales person problem\n",
    "- How to fix dynamic load imbalance\n",
    "- The concept of search overhead\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "480af594",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> Do not forget to execute the cell below before starting this notebook! \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e93809a",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Printf\n",
    "function answer_checker(answer,solution)\n",
    "    if answer == solution\n",
    "        \"🥳 Well done! \"\n",
    "    else\n",
    "        \"It's not correct. Keep trying! 💪\"\n",
    "    end |> println\n",
    "end\n",
    "tsp_check_2(answer) = answer_checker(answer, 2)\n",
    "tsp_check_3(answer) = answer_checker(answer, \"d\")\n",
    "tsp_check_4(answer) = answer_checker(answer, \"a\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1d637cb",
   "metadata": {},
   "source": [
    "## The traveling sales person (TSP) problem\n",
    "\n",
    "\n",
    "### Problem statement\n",
    "\n",
    "Given a graph $G$ with a distance table $C$ and an initial node (i.e. a city) in the graph, compute the shortest route that visits all cities exactly once, without returning to the initial city."
   ]
  },
  {
   "attachments": {
    "fig-tsp-intro.png": {
     "image/png": "[encoded data, 92016 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "bf3dd0e6",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 92016 characters]\"  align=\"left\" width=\"450\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c303dddf",
   "metadata": {},
   "source": [
    "### Sequential algorithm (branch and bound)\n",
    "\n",
    "The sequential algorithm finds a shortest path by traversing the paths tree of the problem. The root of this tree is the initial city. The children of each node in the graph are the neighbour cities that have not been visited on the path so far. When all neighbour cities are already visited, the city becomes a leaf node in the tree.\n",
    "\n",
    "The possile solutions of the problem are the paths from the root of the tree to a leaf node. Note that we assume the children are sorted using the **nearest city first heuristic**. This allows to quickly find a minimum bound for the distance which will be used to prune the remaining paths (see next section). "
   ]
  },
  {
   "attachments": {
    "g25091.png": {
     "image/png": "[encoded data, 78328 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "9797bb92",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 78328 characters]\"  align=\"left\" width=\"525\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da5f5ae",
   "metadata": {},
   "source": [
    "Of course, visiting all paths in the tree is impractical for moderate and large numbers of cities. The number of possible paths might be up to $O(N!)$. Therefore, an essential part of the algorithm is to bound the search by remembering the current minimum distance. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a40faff",
   "metadata": {},
   "source": [
    "### Pruning the search tree\n",
    "\n",
    "The algorithm keeps track of the best solution of all paths visited so far. This allows to skip searching paths that already exceed this value. \n",
    "\n",
    "For example, in the following graph only 3 out of 6 possible routes need to be visted when we cut off the search after the minimum distance is exceeded. (The grey nodes are the ones we don't visit because the minimum distance had been exceeded at the previous node already.)"
   ]
  },
  {
   "attachments": {
    "fig-tsp-pruning.png": {
     "image/png": "[encoded data, 186948 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "ecf89636",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 186948 characters]\"  align=\"left\" width=\"700\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3c78a1d",
   "metadata": {},
   "source": [
    "Note that it is not necessary that the graph be fully connected. Variations of this algorithm work for sparse graphs or directed graphs as well. \n",
    "\n",
    "In the previous example, the shortest route was also the leftmost path in the graph. Although it is more likely that the shortest route be found in the left part of the graph when using the nearest city first heuristic, the solution can be anywhere in the search tree."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85d771de",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b> Example:</b> Look at the following graph and its corresponding search tree. If $x\\leq 15$, the shortest route is the leftmost branch of the search tree. If $x > 16$, the route is situated on the right side of the search tree. \n",
    "</div>"
   ]
  },
  {
   "attachments": {
    "fig-tsp-question-1.png": {
     "image/png": "[encoded data, 209492 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "d1df7878",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 209492 characters]\"  align=\"left\" width=\"800\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fe017a4",
   "metadata": {},
   "source": [
    "## Serial implementation\n",
    "\n",
    "Let's implement the serial algorithm. First, we sort the neighbours according to their distance. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2b70f85",
   "metadata": {},
   "source": [
    "### Nearest-city first heuristic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a50706bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "function sort_neighbors(C)\n",
    "    n = size(C,1)\n",
    "    map(1:n) do i\n",
    "        Ci = C[i,:]\n",
    "        cities = sortperm(Ci)\n",
    "        distances = Ci[cities]\n",
    "        collect(zip(cities,distances))\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2eeecdd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "C = [\n",
    "    0 2 3 2\n",
    "    2 0 4 1\n",
    "    3 4 0 3\n",
    "    2 1 3 0    \n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f769627c",
   "metadata": {},
   "source": [
    "The data structure we will use for the connections table is a matrix of tuples of the form (destination, distance). The tuples are sorted by their distance in ascending order (per start city). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dd0288e",
   "metadata": {},
   "outputs": [],
   "source": [
    "C_sorted = sort_neighbors(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51ed8312",
   "metadata": {},
   "source": [
    "The connections matrix can be indexed by a city. This returns a `Vector{Tuple}}` of all the destinations and their corresponding distances.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00608e1d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "city = 3\n",
    "C_sorted[city]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44025bd5",
   "metadata": {},
   "source": [
    "### Loop over all paths"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c91a99f",
   "metadata": {},
   "source": [
    "Next, we write an algorithm that traverses the whole search tree and prints all the possible paths. The tree is traversed in depth-first order. Before we go to a neighbouring city, we also have to verify that it has not been visited on this path yet. If we reach a leaf node, we print the complete path and continue searching. "
   ]
  },
  {
   "attachments": {
    "fig-tsp-output-1.png": {
     "image/png": "[encoded data, 131352 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "987a0273",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 131352 characters]\"  align=\"left\" width=\"600\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ddc2ec1",
   "metadata": {},
   "outputs": [],
   "source": [
    "function visit_all_paths(C_sorted,city)\n",
    "    # Initialize path\n",
    "    num_cities = length(C_sorted)\n",
    "    path=zeros(Int,num_cities)\n",
    "    hops = 1\n",
    "    path[hops] = city\n",
    "    visit_all_paths_recursive!(C_sorted,hops,path)\n",
    "end\n",
    "function visit_all_paths_recursive!(C_sorted,hops,path)\n",
    "    num_cities = length(C_sorted)\n",
    "    if hops != num_cities\n",
    "        city = path[hops]\n",
    "        next_hops = hops + 1\n",
    "        connections = C_sorted[city]\n",
    "        for (next_city,distance_increment) in connections\n",
    "            # Check if next_city is already on path\n",
    "            already_visited = (next_city in view(path,1:hops))\n",
    "            if ! already_visited\n",
    "                path[next_hops] = next_city\n",
    "                visit_all_paths_recursive!(C_sorted,next_hops,path)\n",
    "            end\n",
    "        end\n",
    "        return nothing\n",
    "    else\n",
    "        println(path)\n",
    "        return nothing\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "723a0f1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "city = 1\n",
    "visit_all_paths(C_sorted,city)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f68093d",
   "metadata": {},
   "source": [
    "### Serial implementation without pruning\n",
    "\n",
    "Now, we add the computation of the minimum distance. At each leaf node, we update the minimum distance. Furthermore, as we add another node to our path, we update the distance of the current path. That makes it necessary to include two more parameters in our recursive algorithm: `distance`, the distance of the current path, and `min_distance`, the best minimum distance found so far. "
   ]
  },
  {
   "attachments": {
    "fig-tsp-output-2.png": {
     "image/png": "[encoded data, 144188 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "881076fa",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 144188 characters]\"  align=\"left\" width=\"600\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5989f0ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "function tsp_serial_no_prune(C_sorted,city)\n",
    "    num_cities = length(C_sorted)\n",
    "    path=zeros(Int,num_cities)\n",
    "    hops = 1\n",
    "    path[hops] = city\n",
    "    distance=0\n",
    "    min_distance = typemax(Int)\n",
    "    tsp_serial_no_prune_recursive!(C_sorted,hops,path,distance,min_distance)\n",
    "end\n",
    "function tsp_serial_no_prune_recursive!(C_sorted,hops,path,distance,min_distance)\n",
    "    num_cities = length(C_sorted)\n",
    "    if hops != num_cities\n",
    "        city = path[hops]\n",
    "        next_hops = hops + 1\n",
    "        connections = C_sorted[city]\n",
    "        for (next_city,distance_increment) in connections\n",
    "            already_visited = (next_city in view(path,1:hops))\n",
    "            if ! already_visited\n",
    "                # Add new city to path\n",
    "                next_distance = distance + distance_increment\n",
    "                path[next_hops] = next_city\n",
    "                min_distance = tsp_serial_no_prune_recursive!(\n",
    "                    C_sorted,next_hops,path,next_distance,min_distance)\n",
    "            end\n",
    "        end\n",
    "        return min_distance\n",
    "    else\n",
    "        # Set new minimum distance in leaf nodes\n",
    "        min_distance = min(distance,min_distance)\n",
    "        #@show path, distance, min_distance\n",
    "        return min_distance\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1be2bfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "city = 1\n",
    "min_distance = tsp_serial_no_prune(C_sorted,city)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb9b8353",
   "metadata": {},
   "source": [
    "### Final serial implementation\n",
    "\n",
    "Finally, we add the pruning to our algorithm. Anytime the current distance exceeds the minimum distance, the search in this path is aborted and continued with another path. "
   ]
  },
  {
   "attachments": {
    "fig-tsp-output-3.png": {
     "image/png": "[encoded data, 135532 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "3c160cf9",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 135532 characters]\"  align=\"left\" width=\"600\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8241e0df",
   "metadata": {},
   "outputs": [],
   "source": [
    "function tsp_serial(C_sorted,city)\n",
    "    num_cities = length(C_sorted)\n",
    "    path=zeros(Int,num_cities)\n",
    "    hops = 1\n",
    "    path[hops] = city\n",
    "    distance=0\n",
    "    min_distance = typemax(Int)\n",
    "    tsp_serial_recursive!(C_sorted,hops,path,distance,min_distance)\n",
    "end\n",
    "function tsp_serial_recursive!(C_sorted,hops,path,distance,min_distance)\n",
    "    # Prune this path if its distance is too high already\n",
    "    if distance >= min_distance\n",
    "        return min_distance\n",
    "    end\n",
    "    num_cities = length(C_sorted)\n",
    "    if hops != num_cities\n",
    "        city = path[hops]\n",
    "        connections = C_sorted[city]\n",
    "        next_hops = hops + 1\n",
    "        for (next_city,distance_increment) in connections\n",
    "            already_visited = (next_city in view(path,1:hops))\n",
    "            if ! already_visited\n",
    "                # Add city to path \n",
    "                next_distance = distance + distance_increment\n",
    "                path[next_hops] = next_city\n",
    "                min_distance = tsp_serial_recursive!(\n",
    "                    C_sorted,next_hops,path,next_distance,min_distance)\n",
    "            end\n",
    "        end\n",
    "        return min_distance\n",
    "    else\n",
    "        # Set new minimum distance in leaf nodes\n",
    "        min_distance = min(distance,min_distance)\n",
    "        #@show path, distance, min_distance\n",
    "        return min_distance\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "998087f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "city = 1\n",
    "min_distance = tsp_serial(C_sorted,city)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f0a5d5c",
   "metadata": {},
   "source": [
    "### Performance\n",
    "Let's have a look at how well our code performs. In the following block, we create another random matrix and time the execution for both versions with and without pruning. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1eb74d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 12 # It is safe to test up to n=12\n",
    "using Random\n",
    "using Test\n",
    "Random.seed!(1)\n",
    "C = rand(1:10,n,n)\n",
    "C_sorted = sort_neighbors(C)\n",
    "city = 1\n",
    "@time min_no_prune = tsp_serial_no_prune(C_sorted,city)\n",
    "@time min_prune = tsp_serial(C_sorted,city)\n",
    "@test min_no_prune == min_prune"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6088ddc9",
   "metadata": {},
   "source": [
    "You can observe that, especially for larger numbers of cities (n=11 or n=12), the performance of the algorithm with pruning is much better than the performance of the algorithm without pruning. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "733ce76e",
   "metadata": {},
   "source": [
    "## Parallel algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6375465",
   "metadata": {},
   "source": [
    "### Where can we extract parallelism ?\n",
    "Unlike the previous algorithms we studied, in this problem we don't know beforehand how much work is performed since we don't know where the pruning cuts off part of the search tree. Still, we want to divide the workload among multiple processes to enhance the performance. "
   ]
  },
  {
   "attachments": {
    "g25091.png": {
     "image/png": "[encoded data, 78328 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "5e853346",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 78328 characters]\"  align=\"left\" width=\"525\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c72df0e",
   "metadata": {},
   "source": [
    "### Option 1\n",
    "\n",
    "The first idea how to parallelize the TSP algorithm is to assign a branch of our search tree to each process. However, as mentioned in an earlier section, the number of branches in the search tree can be up to $O(N!)$. This would require an unfeasibly large amount of proecesses which each do only very little work. "
   ]
  },
  {
   "attachments": {
    "g18142.png": {
     "image/png": "[encoded data, 67600 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "73b7633a",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 67600 characters]\"  align=\"left\" width=\"300\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "052a5676",
   "metadata": {},
   "source": [
    "### Option 2\n",
    "Instead of assigning one branch per worker, we can assign a fixed number of branches to each worker. This way, each worker can perform the pruning within their own subtree and less workers are needed. \n"
   ]
  },
  {
   "attachments": {
    "g18910.png": {
     "image/png": "[encoded data, 58748 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "2701cc6a",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 58748 characters]\"  align=\"left\" width=\"300\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3af7def",
   "metadata": {},
   "source": [
    "### Performance issues\n",
    "\n",
    "#### Load balancing\n",
    "However, this approach has a problem with load balancing. Since we don't know beforehand how much pruning can be done in each subtree, some workers might end up doing less work than others. This uneven distribution of workload leads to some workers being idle, which impairs the speedup. \n",
    "\n",
    "#### Search overhead \n",
    "Another disadvantage of this kind of parallel search is that the pruning is now less effective. The workers each run their own version of the search algorithm and keep track of their local minimum distances. This means that less nodes will be pruned in the parallel version than in the serial version. This is called **search overhead**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9edf5c1",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Question:</b> How many nodes are pruned in total when we assign two branches to each worker? Look at the illustration below.\n",
    "</div>"
   ]
  },
  {
   "attachments": {
    "g20390.png": {
     "image/png": "[encoded data, 118580 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "673fb192",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 118580 characters]\"  align=\"left\" width=\"300\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95e952af",
   "metadata": {},
   "outputs": [],
   "source": [
    "answer = -1 # Replace with any integer value \n",
    "tsp_check_2(answer)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0bc4fdd",
   "metadata": {},
   "source": [
    "In this example, the parallel algorithm prunes less nodes than the serial version because not all workers are able to use the global minimum distance as a pruning bound."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95643cd8",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<b>Question:</b> The previous example described positive search overhead. There is also negative search overhead, resulting in superlinear speedups. Can there be negative search overhead in this parallel TSP algorithm? (Provided the workers communicate the minimum distance with each other)\n",
    "</div>\n",
    "\n",
    "    a) No, because we use the nearest city first heuristic. \n",
    "    b) No, because each worker has to search the whole subtree before the algorithm completes. \n",
    "    c) Yes, because the parallel algorithm does not need to search the whole search tree.\n",
    "    d) Yes, because the global minimum distance can be found more quickly, enabling the parallel version to do more pruning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fff58498",
   "metadata": {},
   "outputs": [],
   "source": [
    "answer = \"x\" # Replace x with a,b,c or d\n",
    "tsp_check_3(answer)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6190b00a",
   "metadata": {},
   "source": [
    "### Option 3: Dynamic load balancing with replicated workers model\n",
    "In our parallel implementation, we will use a coordinator process and several worker processes. The coordinator process (or _master_) searches the tree up to a certain maximum depth _maxhops_. When _maxhops_ is reached, the coordinator creates a job and delegates it to a worker. The workers repeatedly get work from the master and execute it. This is an example of **dynamic load balancing**: the load is distributed among the workers during runtime."
   ]
  },
  {
   "attachments": {
    "fig-tsp-master-worker-2.png": {
     "image/png": "[encoded data, 290908 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "c6d51e00",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 290908 characters]\"  align=\"left\" width=\"600\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "497bad13",
   "metadata": {},
   "source": [
    "This avoids the performance issues of the previous approach: the workers are not idle since they get new work every time they finish a job. To avoid the search overhead, it is also important that the workers communicate the minimum distance each time they update their local minimum distance. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81219a27",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "    <b>Question:</b> To find the right <i> maxhops </i> level is a tradeoff between...\n",
    "</div>\n",
    " \n",
    "    a) Communication overhead (large maxhops) and load imbalance (small maxhops). \n",
    "    b) Search overhead (large maxhops) and load imbalance (small maxhops). \n",
    "    c) the number of workers (large maxhops) and the job size (small maxhops). \n",
    "    d) buffer for the job queue (large maxhops) and idle time of the coordinator process (small maxhops)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cf1ec88",
   "metadata": {},
   "outputs": [],
   "source": [
    "answer=\"x\" #Replace x with a,b,c, or d\n",
    "tsp_check_4(answer)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e4fdbe1",
   "metadata": {},
   "source": [
    "## Implementation of the parallel algorithm \n",
    "\n",
    "First, let's add our worker processes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "354f5abb",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Distributed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f8f6702",
   "metadata": {},
   "outputs": [],
   "source": [
    "if workers() == procs()\n",
    "    addprocs(3)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26188ad3",
   "metadata": {},
   "source": [
    "### Visit paths in parallel\n",
    "\n",
    "Next we will write a simplified version of the algorithm which visits all paths in parallel. The communication between the workers and the coordinator is implemented using a `RemoteChannel`. The coordinator calls the function `generate_jobs` which recursively traverses the search tree until _maxhops_ is reached. Then, it copies the current path to the channel. The workers repeatedly take jobs from the channel until the channel is empty. The workers receive the path from the job and continue the search. When they reach a leaf node in the search tree, they print the path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c0603a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "function visit_all_paths_dist(C_sorted,city,max_hops)\n",
    "    T = typeof((0,Int[]))\n",
    "    jobs_chnl = RemoteChannel(()->Channel{T}(1))\n",
    "    @sync begin\n",
    "        for w in workers()\n",
    "            @spawnat w consume_jobs(C_sorted,jobs_chnl)\n",
    "        end\n",
    "        generate_jobs(C_sorted,city,max_hops,jobs_chnl)\n",
    "        close(jobs_chnl)\n",
    "    end\n",
    "end\n",
    "function generate_jobs(C_sorted,city,max_hops,jobs_chnl)\n",
    "    num_cities = length(C_sorted)\n",
    "    path=zeros(Int,num_cities)\n",
    "    hops = 1\n",
    "    path[hops] = city\n",
    "    generate_jobs_recursive(C_sorted,hops,path,max_hops,jobs_chnl)\n",
    "end\n",
    "function generate_jobs_recursive(C_sorted,hops,path,max_hops,jobs_chnl)\n",
    "    num_cities = length(C_sorted)\n",
    "    if hops == max_hops\n",
    "        @show path\n",
    "        put!(jobs_chnl,(hops,copy(path)))\n",
    "        return nothing\n",
    "    else\n",
    "        city = path[hops]\n",
    "        connections = C_sorted[city]\n",
    "        next_hops = hops + 1\n",
    "        for (next_city,distance_increment) in connections\n",
    "            already_visited = (next_city in view(path,1:hops))\n",
    "            if ! already_visited\n",
    "                path[next_hops] = next_city\n",
    "                generate_jobs_recursive(C_sorted,next_hops,path,max_hops,jobs_chnl)\n",
    "            end\n",
    "        end\n",
    "        return nothing\n",
    "    end\n",
    "end\n",
    "@everywhere function consume_jobs(C_sorted,jobs_chnl)\n",
    "    while true\n",
    "        try\n",
    "            hops,path = take!(jobs_chnl)\n",
    "            consume_jobs_recursive(C_sorted,hops,path)\n",
    "        catch\n",
    "            break\n",
    "        end        \n",
    "    end\n",
    "end\n",
    "@everywhere function consume_jobs_recursive(C_sorted,hops,path)\n",
    "    num_cities = length(C_sorted)\n",
    "    if hops != num_cities\n",
    "        city = path[hops]\n",
    "        connections = C_sorted[city]\n",
    "        next_hops = hops + 1\n",
    "        for (next_city,distance_increment) in connections\n",
    "            already_visited = (next_city in view(path,1:hops))\n",
    "            if ! already_visited\n",
    "                path[next_hops] = next_city\n",
    "                consume_jobs_recursive(C_sorted,next_hops,path)\n",
    "            end\n",
    "        end\n",
    "        return nothing\n",
    "    else\n",
    "        @show path\n",
    "        return nothing\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02801b66",
   "metadata": {},
   "outputs": [],
   "source": [
    "C = [\n",
    "    0 2 3 2\n",
    "    2 0 4 1\n",
    "    3 4 0 3\n",
    "    2 1 3 0    \n",
    "]\n",
    "C_sorted = sort_neighbors(C)\n",
    "city = 1\n",
    "max_hops = 2\n",
    "min_distance = visit_all_paths_dist(C_sorted,city,max_hops)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3bbf609",
   "metadata": {},
   "source": [
    "### How to track the global minimum distance?"
   ]
  },
  {
   "attachments": {
    "fig-tsp-master-worker.png": {
     "image/png": "[encoded data, 290908 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "d85c8e7c",
   "metadata": {},
   "source": [
    "<div>\n",
    "<img src=\"data:image/png;base64,[encoded data, 290908 characters]\"  align=\"left\" width=\"600\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d579e3d",
   "metadata": {},
   "source": [
    "\n",
    "### Simplified example\n",
    "\n",
    "We will demonstrate how the workers communicate the minimum distance with each other with a short example. Each worker generates a random value and updates a globally shared minimum. The variable for the global minimum is stored in a `RemoteChannel`. The buffer size of the channel is 1, such that only one channel can take and put new values to the channel at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0303392",
   "metadata": {},
   "outputs": [],
   "source": [
    "buffer = 1 # Very important\n",
    "min_distance_chnl = RemoteChannel(()->Channel{Int}(buffer))\n",
    "put!(min_distance_chnl,typemax(Int))\n",
    "@sync for w in workers()\n",
    "    @spawnat w begin\n",
    "         sleep(rand(1:3))\n",
    "         min_distance = take!(min_distance_chnl)\n",
    "         @show min_distance\n",
    "         distance = rand(5:10)\n",
    "         @show distance\n",
    "         min_distance = min(distance,min_distance)\n",
    "         put!(min_distance_chnl,min_distance)\n",
    "    end\n",
    "end\n",
    "min_distance = take!(min_distance_chnl)\n",
    "@show min_distance\n",
    "close(min_distance_chnl)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fabefcb",
   "metadata": {},
   "source": [
    "### Final parallel implementation\n",
    "Finally, we add the communication of the minimum distance to the first simplified version. Like in the serial implementation, the workers stop the search if the distance of their current path exceeds the global minimum distance. In total, we need two channels: one to communicate the jobs and one to communicate the minimum distance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96c649b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "function tsp_dist(C_sorted,city,max_hops)\n",
    "    T = typeof((0,Int[],0))\n",
    "    jobs_chnl = RemoteChannel(()->Channel{T}(1))\n",
    "    min_distance_chnl = RemoteChannel(()->Channel{Int}(1))\n",
    "    put!(min_distance_chnl,typemax(Int))\n",
    "    @sync begin\n",
    "        for w in workers()\n",
    "            @spawnat w consume_jobs(C_sorted,jobs_chnl,min_distance_chnl)\n",
    "        end\n",
    "        generate_jobs(C_sorted,city,max_hops,jobs_chnl)\n",
    "        close(jobs_chnl)\n",
    "    end\n",
    "    min_distance = take!(min_distance_chnl)\n",
    "    close(min_distance_chnl)\n",
    "    return min_distance\n",
    "end\n",
    "function generate_jobs(C_sorted,city,max_hops,jobs_chnl)\n",
    "    num_cities = length(C_sorted)\n",
    "    path=zeros(Int,num_cities)\n",
    "    hops = 1\n",
    "    path[hops] = city\n",
    "    distance = 0\n",
    "    generate_jobs_recursive(C_sorted,hops,path,max_hops,jobs_chnl,distance)\n",
    "end\n",
    "function generate_jobs_recursive(C_sorted,hops,path,max_hops,jobs_chnl,distance)\n",
    "    num_cities = length(C_sorted)\n",
    "    if hops == max_hops\n",
    "        #@show path, distance\n",
    "        put!(jobs_chnl,(hops,copy(path),distance))\n",
    "        return nothing\n",
    "    else\n",
    "        city = path[hops]\n",
    "        connections = C_sorted[city]\n",
    "        next_hops = hops + 1\n",
    "        for (next_city,distance_increment) in connections\n",
    "            already_visited = (next_city in view(path,1:hops))\n",
    "            next_distance = distance + distance_increment\n",
    "            if ! already_visited\n",
    "                path[next_hops] = next_city\n",
    "                generate_jobs_recursive(\n",
    "                        C_sorted,next_hops,path,max_hops,jobs_chnl,next_distance)\n",
    "            end\n",
    "        end\n",
    "        return nothing\n",
    "    end\n",
    "end\n",
    "@everywhere function consume_jobs(C_sorted,jobs_chnl,min_distance_chnl)\n",
    "    while true\n",
    "        min_distance = take!(min_distance_chnl)\n",
    "        put!(min_distance_chnl,min_distance)\n",
    "        try\n",
    "            hops,path,distance = take!(jobs_chnl)\n",
    "            consume_jobs_recursive(C_sorted,hops,path,min_distance_chnl,distance,min_distance)\n",
    "        catch\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "end\n",
    "@everywhere function consume_jobs_recursive(C_sorted,hops,path,min_distance_chnl,distance,min_distance)\n",
    "    if distance >= min_distance\n",
    "        return min_distance\n",
    "    end\n",
    "    num_cities = length(C_sorted)\n",
    "    if hops != num_cities\n",
    "        city = path[hops]\n",
    "        connections = C_sorted[city]\n",
    "        next_hops = hops + 1\n",
    "        for (next_city,distance_increment) in connections\n",
    "            already_visited = (next_city in view(path,1:hops))\n",
    "            next_distance = distance + distance_increment\n",
    "            if ! already_visited\n",
    "                path[next_hops] = next_city\n",
    "                min_distance = consume_jobs_recursive(\n",
    "                        C_sorted,next_hops,path,min_distance_chnl,next_distance,min_distance)\n",
    "            end\n",
    "        end\n",
    "        return min_distance\n",
    "    else\n",
    "        min_distance = take!(min_distance_chnl)\n",
    "        min_distance = min(min_distance,distance)\n",
    "        put!(min_distance_chnl,min_distance)\n",
    "        #@show path, distance, min_distance\n",
    "        return min_distance\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67888155",
   "metadata": {},
   "outputs": [],
   "source": [
    "city = 1\n",
    "max_hops = 2\n",
    "min_distance = tsp_dist(C_sorted,city,max_hops)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48b54da1",
   "metadata": {},
   "source": [
    "## Testing the parallel implementation\n",
    "Next, we will test the correctness and performance of our parallel implementation by comparing the results of the parallel algorithm to the results of the serial algorithm for multiple problem instances. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d950905",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 18 # Safe to run up to 18\n",
    "using Random\n",
    "Random.seed!(1)\n",
    "C = rand(1:10,n,n)\n",
    "C_sorted = sort_neighbors(C)\n",
    "city = 1\n",
    "T1 = @elapsed min_serial = tsp_serial(C_sorted,city)\n",
    "max_hops = 2\n",
    "P = nworkers()\n",
    "TP = @elapsed min_dist = tsp_dist(C_sorted,city,max_hops)\n",
    "println(\"Sequential time = \", T1)\n",
    "println(\"Parallel time = \", TP)\n",
    "println(\"Speedup = \", T1/TP)\n",
    "println(\"Optimal speedup = \", P)\n",
    "println(\"Efficiency = \", 100*(T1/TP)/P, \"%\")\n",
    "@test min_serial == min_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c789dc7a",
   "metadata": {},
   "source": [
    "# License\n",
    "\n",
    "\n",
    "\n",
    "This notebook is part of the course [Programming Large Scale Parallel Systems](https://www.francescverdugo.com/XM_40017) at Vrije Universiteit Amsterdam and may be used under a [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) license."
   ]
  }
 ],
 "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
}
