Djikstra_Bellman_Ford_Heapq.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Well-known graph traversal implimentations (simplified):\n",
    "\n",
    "#### 1- Djikstra Algorithm to find the minimum shortest path. O(V**2)\n",
    "#### 2- Priority Heap \n",
    "#### 3- Djikstra Algorithm improved with the priority heap so that the running time is O(V*logE)  \n",
    "#### 4- Bellman Ford Algorithm to find the minimum shortest path. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  inf 71.52 60.28 54.49 42.37 64.59 43.76 89.18 96.37 38.34]\n",
      " [79.17   inf 56.8  92.56  7.1   8.71  2.02 83.26 77.82 87.  ]\n",
      " [97.86 79.92   inf 78.05 11.83 63.99 14.34 94.47 52.18 41.47]\n",
      " [26.46 77.42 45.62   inf  1.88 61.76 61.21 61.69 94.37 68.18]\n",
      " [35.95 43.7  69.76  6.02   inf 67.06 21.04 12.89 31.54 36.37]\n",
      " [57.02 43.86 98.84 10.2  20.89   inf 65.31 25.33 46.63 24.44]\n",
      " [15.9  11.04 65.63 13.82 19.66 36.87   inf  9.71 83.79  9.61]\n",
      " [97.65 46.87 97.68 60.48 73.93  3.92 28.28   inf 29.61 11.87]\n",
      " [31.8  41.43  6.41 69.25 56.66 26.54 52.32  9.39   inf 92.93]\n",
      " [31.86 66.74 13.18 71.63 28.94 18.32 58.65  2.01 82.89   inf]]\n"
     ]
    }
   ],
   "source": [
    "# prepare a densely connected test matrix with random numbers\n",
    "# to represent a densely connected directed cyclic graph, initialize the distances as 0 from each vertex to itself\n",
    "\n",
    "import numpy as np\n",
    "np.random.seed(0)\n",
    "def test_graph(n): # n vertices\n",
    "    test_graph_set = np.random.rand(n,n)*100\n",
    "    # 0 the source so that the node is acyclic to itself \n",
    "    for i in range(n):\n",
    "        test_graph_set[i][i] = np.inf\n",
    "        # test_graph_set[i][:i] = np.inf  # unmark if a tree structure is needed\n",
    "    test_graph_set = np.around(test_graph_set, decimals = 2)\n",
    "    return test_graph_set\n",
    "\n",
    "test_set = test_graph(10) \n",
    "print(test_set[:10][:10])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Naive Djikstra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.   54.8  51.52 48.39 42.37 44.27 43.76 40.35 69.96 38.34]\n",
      "95.7 µs ± 6.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# write a naive DIJKSTRA's algorithm to find the shortest distance from a source vertex to all other vertices (O(v**2))\n",
    "def graph_traverse_DJ_path(graph, source): \n",
    "    distance_list = [np.inf] * len(graph)\n",
    "    visited_vertices_list = [False] * len(graph)\n",
    "    distance_list[source] = 0\n",
    "    for each_vertex in range(len(graph)):\n",
    "        min_value = np.inf\n",
    "        for i in range(len(distance_list)):\n",
    "            if distance_list[i] < min_value and visited_vertices_list[i] == False:\n",
    "                min_value = distance_list[i]\n",
    "                min_at_index = i\n",
    "        visited_vertices_list[min_at_index] = True\n",
    "        for i in range(len(distance_list)):\n",
    "            if visited_vertices_list[i] == False and graph[min_at_index][i] > 0 and \\\n",
    "            distance_list[min_at_index] + graph[min_at_index][i] < distance_list[i]: \n",
    "                distance_list[i] = distance_list[min_at_index] + graph[min_at_index][i]\n",
    "    return np.around(distance_list, decimals = 2)\n",
    "\n",
    "print(graph_traverse_DJ_path(test_set, 0))\n",
    "%timeit graph_traverse_DJ_path(test_set, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parent index: 0 left child index: 1 right child index: 2\n",
      "parent index: 1 left child index: 3 right child index: 4\n",
      "parent index: 2 left child index: 5 right child index: 6\n",
      "parent index: 3 left child index: 7 right child index: 8\n",
      "parent index: 4 left child index: 9 right child index: 10\n"
     ]
    }
   ],
   "source": [
    "# visualize the parent - child nodes' index relationship in a heap:\n",
    "nodes_indices = np.arange(0,5)\n",
    "index = 0\n",
    "right_child = 0\n",
    "while index < len(nodes_indices):\n",
    "    left_child = index *2 + 1\n",
    "    right_child = index *2 + 2\n",
    "    print('parent index:', index, 'left child index:', left_child, 'right child index:', right_child)\n",
    "    index += 1\n",
    "\n",
    "# percolate up: while index i > 0: move from last index to top and switch the parent and left child, update i.\n",
    "# percolate down: while index i < len//2: find min child of that node, switch the parent and mc, update i as mc "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Priority Heap (minimum at first index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 8, 3, 27, 8, 6, 5245, 245, 78, 26, 725, 24, 9]\n",
      "[[0.0, 'hello'], [10.02, 'first'], [10.32, '0'], [224.09, 'dog'], [186.76, 'oink'], [97.87, 'bear']]\n",
      "0.0 <class 'float'> hello <class 'str'>\n"
     ]
    }
   ],
   "source": [
    "# create a heapq to later use it in the djikstra's algorithm\n",
    "class Heapq:\n",
    "    def __init__(self):\n",
    "        self.heap = []\n",
    "    \n",
    "            \n",
    "    def push(self,k): # append k to the array and percolate up with the last index\n",
    "        self.heap.append(k)\n",
    "        index_of_last_item = len(self.heap)-1 \n",
    "        self.percUp(index_of_last_item) # last item's index\n",
    "\n",
    "             \n",
    "    def percUp(self,i): # commonly known as \"siftdown method\"\n",
    "        while i > 0: # smallest index need to be 1 since first index is 0\n",
    "            if self.heap[i] < self.heap[(i-1) // 2]: # if child is smaller than its parent, swap\n",
    "                tmp = self.heap[(i-1) // 2]\n",
    "                self.heap[(i-1) // 2] = self.heap[i]\n",
    "                self.heap[i] = tmp\n",
    "            i = (i-1) // 2\n",
    "\n",
    "            \n",
    "    def pop(self):\n",
    "        if len(self.heap) == 0: return None\n",
    "        if len(self.heap) == 1: \n",
    "            tmp = self.heap.pop()\n",
    "            return tmp\n",
    "        retval = self.heap[0] # first item\n",
    "        last_item = self.heap.pop() # delete last item and assign it to a variable\n",
    "        self.heap[0] = last_item  # swap first item with last item\n",
    "        self.percDown(0)\n",
    "        return retval\n",
    "\n",
    "\n",
    "    def percDown(self,i): # commonly known as \"siftup method\"\n",
    "        # leftchild index is 2*i + 1, rightchild index is 2*i + 2, \n",
    "        # rightchild might not exist\n",
    "        while i < (len(self.heap) // 2):\n",
    "            mc = self.minChild(i)\n",
    "            if self.heap[i] > self.heap[mc]:\n",
    "                tmp = self.heap[i]\n",
    "                self.heap[i] = self.heap[mc]\n",
    "                self.heap[mc] = tmp\n",
    "            i = mc   \n",
    "            \n",
    "            \n",
    "    def minChild(self,i):\n",
    "        left_child_index = (i*2) + 1\n",
    "        right_child_index = left_child_index + 1\n",
    "        # if there is no right child, return the left child's index\n",
    "        if (left_child_index + 1) == len(self.heap): return left_child_index\n",
    "        elif self.heap[left_child_index] < self.heap[right_child_index]: # choose the smaller child \n",
    "            return left_child_index\n",
    "        else:\n",
    "            return right_child_index\n",
    "\n",
    "\n",
    "    def heapify(self,alist): \n",
    "        self.heap = alist[:]\n",
    "        n = len(self.heap)\n",
    "        for i in reversed(range(n//2)): # 4, 3, 2, 1, 0 and 0 is also included.\n",
    "            self.percDown(i)\n",
    "        return self.heap\n",
    "    \n",
    "    \n",
    "         \n",
    "if __name__== \"__main__\":\n",
    "            \n",
    "    bh_test = Heapq()\n",
    "    test_set2 = [24,26,6,27,725,3,5245,245,78,8,8,0,9] \n",
    "    bh_test.heapify(test_set2)\n",
    "    print(bh_test.heap)\n",
    "    test_set3 = [[0.0, 'hello'], [10.02, 'first'], [97.87, 'bear'], [224.09, 'dog'], [186.76, 'oink'], [10.32, '0']]\n",
    "    bh_test.heapify(test_set3)\n",
    "    print(bh_test.heap)\n",
    "    value, text = bh_test.pop()\n",
    "    print(value, type(value), text, type(text))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the heap with the built-in heap and a python built-in quick sort to see its efficiency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "546 µs ± 41.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "[5, 10, 14, 12, 12, 17, 26, 25, 24, 34, 30, 20, 48, 111, 42, 29, 33, 149, 43, 45]\n",
      "29.3 µs ± 1.91 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "[5, 10, 14, 12, 12, 17, 26, 25, 24, 34, 30, 20, 48, 111, 42, 29, 33, 149, 43, 45]\n",
      "125 µs ± 2.92 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "[ 1  1  4  4  5  8 13 14 18 21 22 23 25 27 29 30 33 36 43 47]\n"
     ]
    }
   ],
   "source": [
    "# check the efficiency of the heap above and compare your heapq above with a python's built in module\n",
    "import heapq # python built-in heap \n",
    "import numpy as np\n",
    "np.random.seed(0)\n",
    "testing_list = list(np.random.randint(1,1000, size = 500))\n",
    "\n",
    "# my_heap object below takes the list and heapifies it and outputs the object as priority heap:\n",
    "my_heap = Heapq()\n",
    "%timeit my_heap.heapify(testing_list)\n",
    "print(my_heap.heap[:20])\n",
    "\n",
    "# using python built-in heap.heapify to convert a list into priority heap in place\n",
    "python_heap_sorted = heapq.heapify(testing_list) \n",
    "%timeit heapq.heapify(testing_list) \n",
    "print(testing_list[:20])\n",
    "\n",
    "# Additional: The python built-in sorting nlogn algorithm to sort the whole list\n",
    "quick_sorted = np.sort(list(np.random.randint(1,1000, size = 500)))\n",
    "%timeit np.sort(list(np.random.randint(1,1000, size = 500)))\n",
    "print(quick_sorted[:20])\n",
    "\n",
    "# Conclusion: The heapq class written above is working very slowly. Python built-in heapq class is the fastest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Improved Djikstra with Heapq "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{0: 26.46, 1: 33.959999999999994, 2: 45.62, 3: 0, 4: 1.88, 5: 18.689999999999998, 6: 22.919999999999998, 7: 14.77, 8: 33.42, 9: 26.64}\n",
      "126 µs ± 4.02 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# Keep the distances in a dictionary where the keys are the name of vertices and the values are the total \n",
    "# minimum distances traversed from source node to other nodes\n",
    "\n",
    "import numpy as np\n",
    "import heapq\n",
    "\n",
    "def djikstra_with_heapq(graph, starting_vertex):\n",
    "    distances = {vertex: np.inf for vertex in range(len(graph))} # {'0': inf, '1': inf, '2': inf,...}\n",
    "    distances[starting_vertex] = 0 # {'0': 0, '1': inf, '2': inf,...} it is made acyclic\n",
    "    visited_vertices_dict = [False] * len(graph) # {'0': 'False', '1': 'False',..} to make sure each vertex is visited once\n",
    "    # loop over vertices\n",
    "    for k in distances.keys(): # for each vertex in hash_table, append the connection list to heapq\n",
    "        # find the min distance by heapifying the distances list, \n",
    "        # the values and keys need to switch places since values are going to be compared and sorted\n",
    "        distances_as_heap = [[value, key] for key, value in distances.items()]\n",
    "        heapq.heapify(distances_as_heap) # heap: [[0.0, '0'], [10.32, '8'], [95.01, '6'],...]\n",
    "        # pop the first item that is not marked as visited\n",
    "        while len(distances_as_heap) > 0: # '0' '1' '2' ...\n",
    "            min_distance, index_at_min_dist = distances_as_heap.pop(0) \n",
    "            if visited_vertices_dict[index_at_min_dist] == False: break\n",
    "        visited_vertices_dict[index_at_min_dist] = True\n",
    "        # loop over the hash table keys and check if the registered distances[vertex] has a shorter connection\n",
    "        for i in range(len(graph)):\n",
    "            if visited_vertices_dict[i] == False and graph[index_at_min_dist][i] > 0 and \\\n",
    "            distances[index_at_min_dist] + graph[index_at_min_dist][i] < distances[i]: \n",
    "                distances[i] = distances[index_at_min_dist] + graph[index_at_min_dist][i]          \n",
    "    return distances\n",
    "print(djikstra_with_heapq(test_set, 3))\n",
    "%timeit djikstra_with_heapq(test_set, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adjacency list representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0, 5, 64.59], [0, 6, 43.76], [0, 7, 89.18], [0, 8, 96.37], [0, 9, 38.34], [1, 0, 79.17], [1, 1, inf], [1, 2, 56.8], [1, 3, 92.56], [1, 4, 7.1]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "# turn the matrix into an adjacency list\n",
    "adj_list = [[i, u, test_set[i][u]] for i in range(len(test_set)) for u in range(len(test_set))] \n",
    "# view the adjacency list\n",
    "print(adj_list[5:15])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bellman-Ford Algorithm for Graph Traversal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[26.46, 33.959999999999994, 39.82, 0, 1.88, 18.689999999999998, 22.919999999999998, 14.77, 33.42, 26.64]\n",
      "675 µs ± 38 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "120 µs ± 3.95 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "134 µs ± 3.89 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# Bellman-Ford Algorithm to traverse the test graph\n",
    "# print(test_set) to view the matrix representation if needed\n",
    "\n",
    "def graph_bellman_ford(adj_list, source): \n",
    "    # Initialize distances list from source vertex to other vertices\n",
    "    dist = [np.inf] * len(test_set)\n",
    "    # Make it acyclic to itself\n",
    "    dist[source] = 0 \n",
    "    # Loop over V-1 times and update the distance list where V is the number of vertices of the test_set (matrix)\n",
    "    for i in range(len(test_set) - 1): \n",
    "        for u, v, w in adj_list: \n",
    "            if dist[u] != np.inf and w != 0 and dist[u] + w < dist[v]: \n",
    "                dist[v] = dist[u] + w \n",
    "\n",
    "    # Check for negative-weight cycles going over it one more time to see if there is any update.  \n",
    "    for u, v, w in adj_list: \n",
    "        if dist[u] != np.inf and w != 0 and dist[u] + w < dist[v]: \n",
    "            print (\"Graph contains negative weight cycle\")\n",
    "            return\n",
    "    \n",
    "    return dist\n",
    "        \n",
    "print(graph_bellman_ford(adj_list, 3)) \n",
    "# compare naive djikstra, improved djikstra, and bellman-ford\n",
    "%timeit graph_bellman_ford(adj_list, 3) \n",
    "%timeit graph_traverse_DJ_path(test_set, 3)\n",
    "# since toy matrix is very small the performance of the heap is not visiable\n",
    "%timeit djikstra_with_heapq(test_set, 3) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes:\n",
    "Whether the graph is directed, cyclic, acyclic, or undirected does not matter. Dijkstra's algorithm would work. \n",
    "\n",
    "Dijkstra’s algorithm doesn’t work for graphs with negative weight edges. \n",
    "\n",
    "For graphs with negative weight edges, Bellman–Ford algorithm can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}