Functional and Logic Programming - 3rd Year 1st Semester - 1st Project (TP1)

Group T12G02

Guilherme Moura Pinho Andrade Silva, up202205298@up.pt

Valentina Pereira Cadime, up202206262@up.pt

Contribution

Guilherme Silva (50%)

• MinHeap implementation of the priority queue data structure
• cities function and its auxiliary rmDup
• areAdjacent function
• distance function
• adjacent function
• pathDistance function and its auxiliary accPath
• rome function
• isStronglyConnected function and its auxiliary dfs
• shortestPath function and its auxiliaries dijkstra, lookupIndex, fromMaybe, dijkstraLoop, edgeRelax, reconstructPaths, indexToCity and reconstruct
• Documentation of MinHeap, cities, rmDup areAdjacent, distance, adjacent, pathDistance, accPath, shortestPath, dijkstra, lookupIndex, fromMaybe, dijkstraLoop, edgeRelax, reconstructPaths, indexToCity and reconstruct in TP1
• shortestPath section in README

Valentina Cadime (50%)

• travelSales function and its auxiliaries graphToMatrix, fillTable, addMaybe, findPathByIndex and getPath
• Documentation of rome, isStronglyConnected, dfs, travelSales, graphToMatrix, fillTable, addMaybe, findPathByIndex and getPath in TP1
• travelSales section in README

Implementation of the shortestPath function

This function's main goal is to compute all shortest paths connecting two given cities in a given graph.

1) Preparation

MinHeap implementation of the priority queue data structure, as per recommended, to use in Dijkstra's algorithm

• Functions: emptyMinHeap, insert, heapifyUp, extractMin, heapifyDown, swap

• Initial list of tuples-based MinHeap revealed limitations in the element-swapping function - swap - in terms of complexity, resulting in O(n), causing insertion and extraction operations to be O(n log n), expectedly providing the same complexity as using a simple list of tuples to sort vertices by smallest distance. Therefore, we decided to use an Array-based MinHeap, which guarantees O(1) complexity for the swap function, and consequently O(log n) for insertion and extraction operations.

data MinHeap a = MinHeap (Data.Array.Array Int a) Int Int

• Additionally, even with the Array-based MinHeap, due to the Array's immutability, the insert function had O(n log n) complexity, due to the need to create a new array, copy the old array into it, and add the new element - O(n). To solve this, in a similar approach to how dynamic arrays like std::vector in C++ handle resizing, we double the capacity of the array only when the current size reaches its limit. For that we added 2 integer values to the declared MinHeap data-type - currentSize and capacity (as shown above), and changed the insert function's implementation:

insert :: (Ord a) => a -> MinHeap a -> MinHeap a
insert x (MinHeap arr currSize capacity)
  | currSize == capacity = insert x (MinHeap newArr currSize (2*capacity))
  | otherwise = MinHeap (heapifyUp currSize (arr Data.Array.// [(currSize, x)])) (currSize+1) capacity
  where
    -- Double the capacity of the array when it's full (to avoid resizing the array for each insertion). Fill empty space with 'undefined'
    newArr = Data.Array.listArray (0, 2*capacity-1) (Data.Array.elems arr ++ replicate capacity undefined)

• This way, we amortize the insert's function complexity, giving it an average time complexity of O(1), as the Array is resized only when it's full (less frequently), and the worst-case time complexity is O(n log n), where n is the number of elements in the heap.

• The heapifyUp function is called after inserting an element into the MinHeap - by the insert function - to maintain the MinHeap property, i.e. the parent of the inserted element is smaller than the element itself. Therefore, the function compares the element with its parent and swaps them if the parent is greater than the element, recursively calling itself until the MinHeap property is satisfied.

• The extractMin function, if it isn't empty (in which case it returns Nothing), returns a tuple with the root of the MinHeap (the city with the smallest distance to the source, in the context of Dijkstra's algorithm) and the resulting MinHeap after extracting that element, and reorganizes the heap to maintain the MinHeap property - heapifyDown.

• The heapifyDown function compares the new root its children, and swaps it with the smallest child if it's greater than it, recursively calling itself until all nodes are smaller than both children (the MinHeap property is satisfied).

• The swap function swaps two elements in the MinHeap array.

• The emptyMinHeap function initializes an empty MinHeap with a capacity of 1 and a current size of 0.

Transformation of the graph into an adjacency list

• The graph (Roadmap) is represented as a list of edges, where each edge is a tuple of two cities and the distance between them. To cache the adjacent cities of each city in dijkstra, we transform the graph into an adjacency list, as recommended by the guidelines, by calling the adjacent function on each city in the graph (listed with the cities function). Each city is associated with a list of tuples, each containing the adjacent city and the distance between them. To enhance readability, we also added a type alias AdjList for the adjacency list (as shown in the guidelines):

type AdjList = [(City,[(City,Distance)])]
-- ...
-- In the 'dijkstra' function
adjList = [(city, adjacent graph city) | city <- cities graph]

2) Initialization Phase

As per the guidelines, if the source city is the same as the destination city, the shortest path is the city itself. Therefore, we added a guard clause to the shortestPath function to handle this case. Otherwise, it calls the dijkstra function to compute the shortest paths from the given source city to the given destination city.

To code the shortestPath function, we tried to follow the Dijkstra's algorithm pseudocode as closely as possible. So, we start by initializing the distance array, the predecessors array, and the MinHeap:

• The distance array stores the shortest distance from the source city to each city in the graph. In the initial state, the source city is associated with a distance of 0, and all other cities are associated with an infinite distance.
• The predecessors array stores the predecessor cities of each city in the shortest path from the source city. In the initial state, every city is associated with an empty list.
• The MinHeap stores the cities in the graph sorted by their distance from the source city. In the initial state, the source city is inserted with a distance of 0.

initPQ = insert (0, srcIdx) emptyMinHeap
initDist = Data.Array.listArray (0, length allCities-1) (repeat maxBound) Data.Array.// [(srcIdx, 0)]
initPred = Data.Array.listArray (0, length allCities-1) (repeat [])

It is important to note the order of the tuples stored in the priority queue. Opposite to how they are stored in the distance array, the distances are the first element, because the sorting done by the MinHeap is based on the Ord instance of the tuple, which compares the first element of the tuple, and then the second element if the first elements are equal. The distance being the first element naturally leads to the correct ordering of the priority queue in Dijkstra's algorithm, as well as avoids the creation of a custom comparison function for the MinHeap or a custom Ord instance for the tuple.

Until the very last step of building the paths, we use the cities' indexes to represent them in the data structures throughout the algorithm. To convert the cities' names to their respective indexes, we use the lookupIndex function, which returns the index of a city in the list of all cities.

3) Loop Phase

After the initialization, we call the dijkstraLoop function to compute the shortest paths from the source city to the destination city and get the updated distance and predecessor arrays.

In the dijkstraLoop function, we extract the city with the smallest distance (root) from the MinHeap and relax its adjacent cities (edgeRelax function), which updates the priority queue, the distances array and the predecessors array. Then we keep recursively calling dijkstraLoop with those updated values until the MinHeap is empty.

The fromMaybe function is used to extract the value from a Maybe type, returning a default value (empty list) if it's Nothing. This is used to extract the adjacent cities of the current city from the adjList, as the lookup function returns a Maybe type.

The edge relaxation process - in the edgeRelax function - consists of comparing the current city with each of its neighbors. More concretely, we compare the distance from the source city to the current city plus distance from the current city to the adjacent city with the distance from the source city to the adjacent city.

• If the former is smaller, we update the distance array with the new distance, and reset the neighbor's predecessor array with the current city as its only predecessor. We also insert the adjacent city into the MinHeap with the new distance.
• If the values are equal, we append the current city to the list of predecessors of the adjacent city. As the neighbor already has a distance, it means that it was already inserted into the MinHeap, so we don't need to insert it again.
• If the former is greater, we don't need to update the distance array nor the predecessor array.

4) Termination Phase

This phase is reached when the MinHeap is empty, meaning that all cities have been relaxed. The dijkstra function then calls the reconstructPaths function to build the shortest paths from the source city to the destination city using the predecessors array. This function, given the paths containing the indexes of cities from the destination city to the source city - calculated the by reconstruct function - reverses them, converts them into their names and returns them as paths from the source city to the destination city.

In the reconstruct function, the paths are built from the destination city to the source city. Initially it's called with the destination city. So, we start by adding the that city to the path and recursively calling reconstruct with the predecessors of the destination city, and then the city after that, and so on, until we reach the source city (base case). The function returns a list of paths made up of indexes, where each path is a list of city indexes from the destination city to the source city.

5) Conclusion

One characteristic of this Dijkstra's algorithm implementation in Haskell is the prevalent use of lists (for example, in the definition of the Adjacency List, as per the guidelines), which increases the time complexity of the algorithm - when compared with a class implementation in C++, where the adjacent vertices can be a member variable, with constant lookup time. This is because the time complexity of the lookup and !! function is O(n), where n is the length of the list. So even though the adjList variable exists, and the adjacent doesn't need to be repeatedly called, the lookups are still time consuming. This lookup is used both for finding the adjacent cities of a city - in dijkstraLoop - and for finding the index of a city in the list of all cities - as in reconstructPaths. This could be improved by using a more efficient data structure, like a Map, but would require a significant change in the return values of other functions, new imports, and other changes that are not allowed in this project.

Given the process of finding the shortest path between the two cities, while filling the predecessors array, and then reconstructing the path backwards, the most time-consuming step is the latter - the reconstructPaths function - which not only reconstructs all the paths but also converts the resulting list of index paths to paths with city names. The time complexity of this function is O(P * V^2), where P is in the number of shortest paths, V is the number of cities in the graph. In the worst case, P is V! - in a complete graph - which is the number of permutations of V cities. Therefore this is the complexity of the shortestPath function.

Implementation of the travelSales function - Traveling Salesperson Problem

This function's main goal is to get the path that connects all the nodes of a graph, starting and ending at the same place. Though the cost of travel must meet the minimum cost to traverse the nodes.

As mentioned in the guidelines, it isn't mandatory to showcase all paths possible, only one is enough to prove that the function achieves its goal. Therefore, the implementation of this function will be centered on starting and ending the path on the first node obtained from the graph variable, when we traverse the entire graph. We will call the first node 0 (from its index in the Arrays), and the rest of the nodes will consequently follow that numeration.

Ex.: 0 -> 1 -> 3 -> 2 -> 0 or [0,1,3,2,0]

1) The Basics

This function is required to have a dynamic programming approach for the optimal solution of the Traveling Salesperson Problem. Thus, certain auxiliary functions and 'attributes' are needed, such as:

• 1º The dp table with the minimum cost to visit x,y,z,... cities (mask) ending in city 0, depending in which city (index) we currently are.

• 2º The fillTable function to fill that table depending on the graph's edges and their weights.

• 3º The findPathByIndex function that uses the filled table to obtain the minimum cost path from the source, passing through all cities and returning to the source. The method used is called backtracking, where we converge to multiple base cases and then choose the more viable one. After that, we simply follow the path that led to that base case backwards.

• 4º An additional function, graphToMatrix, to transform the original graph (list of edges) into an 2-dimensional array where each pair of indexes (start,end) indicates the direct distance between cities - the weight of the edge connecting cities start and end. That matrix is stored in adjMatrix.

• 5º A list of cities to obtain the correct name of the city depending on the index given. We will reuse the function cities for that matter.

Besides, throughout the implementation, we will use helper function addMaybe to calculate the addition of Maybe Distance, the data type to describe the distance between two cities. This will make the program more readable and modular.

2) The Beginning - Filling the dp table with the minimum cost

Firstly, we have to consider what is a valid graph, and what the function should return in turn.

• If a graph has less or equal than 2 cities, and knowing that the salesperson can't use the same edge twice, then the salesperson can't go back to the original city, hence returning immediately [].
• Otherwise, the function may proceed into filling the dp table.

The dp table consists of a 2D array: with the mask and the index, as described above. We use the mask to include all the cities that we still need to visit before city 0 - which saves space when compared with using a list for the same purpose. Each bit in the mask represents a city's index.

For example, the value of dp(12,1) represents the minimum cost to visit from city 1 to city 0 passing through cities 2 and 3.

Notes: 12 == 1100 in binary, which means that cities 3 and 2 (since bits 3 and 2 are set) still need to be visited. In this case, we don't set the rightmost bit (the bit representing that we need to visit city 0), because it is unnecessary, since we are immediately adding the distance from the final city in the path to city 0, when only that city is left to explore - base case.

We fill the dp table from the simplest case (i.e. one city to another) to where all cities have to be visited. The result is the sum of the minimum cost from visiting all the other cities - except one of them - plus the distance from that one city to city 0. Of course, choosing the best last city to visit is very difficult, since any difference in path can give different and even impossible (Nothing) distances. That's why we compare all the distances given from each city that isn't the current city, and choose the minimal (we exclude the Nothing distances if there is any Just distance, and among the Just distances, the smallest one is chosen).

3) Choosing the correct path

After the dp table is completed, all that it's left to do is follow the minimal path. We know that we have to start from the highest point and reach the base case recursively. That highest point being all cities must be visited and we have to start at 0 (reaching 0 is assumed from the calculations of the table).

Consequently, getPath starts finding the path by having the mask all set (1 in every bit) - (1 'Data.Bits.shiftL' n) - 1 - and the start city bit be 0. In more detail, getPath checks for all cities except 0 (the current index), which one has the smallest distance possible if we had to visit that city right now.

We exclude all the Nothing paths, since that value indicates that it will be impossible to reach 0. Besides, we store the final distance in a tuple with its respective city index, just so we can save some trouble of calculating the distances again when finding the smallest one (memoization).

The relevant code is:

allCitiesToVisit = [(dist, j) | j <- [0..(n-1)], Data.Bits.testBit mask j, j /= indexNow, let dist = addMaybe (adjMatrix Data.Array.! (indexNow, j)) (dp Data.Array.! (newMask, j)), dist /= Nothing]

The relative distance from the current city to 0 is already calculated on the dp table, and the distance between the current index and the city we would like to visit is also accessible through the adjMatrix. So allCitiesToVisit will be easily filled with all the valid candidates to be the city to visit next.

Finally, Data.List.minimumBy decides which city to visit next and that city is added on the final list. We still need to know the remaining movements the salesperson must make to reach city 0, so we use recursion to get one step closer to the base case - having only one city left. The recursion is simply taking the chosen city out of the mask and putting it as the next index, to find the next city to visit, and so on.

When it finally reaches the '1 remaining city' which is equal to the index, then we can add that element to the list alongside the last city to travel - city 0. In the end, the recursion will return the list with all the elements in order.

Remember, that if there's no possible path, in other words, allCitiesToVisit is null, then getPath will return an empty list without the 0 (initial city), since it can't reach base case.

4) The End

Ultimately, all that travelSales has to do is map the name of the cities with their respective index, and return the list.

The time complexity of the travelSales function is O(2^n * n^2), where n is the number of cities in the graph. This is because the fillTable function has a time complexity of O(2^n * n^2), and the findPathByIndex function has a time complexity of O(n^2), where n is the number of cities in the graph. Therefore, the filling the dp table is the most time consuming step.