Table des matières

# Important instructions

**Do not change the signatures of the requested functions**(function names, function parameters...): otherwise you risk losing points.- You must
**clearly answer the questions asked**in your python file as**comments**. - Your file to be submitted must have the name
**votreNom_votrePrenom_votreNumEtu.py**(your student number starts with 23 with 7 digits). - You must submit your file to
**your own group**on Edunao - We will evaluate your work done by yourself and in your group's classroom!

Please respect all these instructions, otherwise your work will not be evaluated (or not evaluated correctly).

# Goals of Lab session 2

Shortest path algorithms are a family of algorithms designed to solve the shortest path problem. The shortest path problem can be defined for undirected or directed graphs. Shortest path algorithms have many applications in our real world, such as GPS systems, Internet routing protocols, power distribution networks, enemy AI navigation in video games, contact tracing and disease spread prediction, etc.

In lectures and tutorials, you have seen that the shortest path / SP algorithm (Dijkstra) is used to find the shortest paths between a given vertex and all others in a graph with only positive weights. The Bellman-Ford algorithm can work for negative weights as well as for the detection of absorbing circuits, but with higher complexity. In this TP session (3h), we'll look at a new problem: we want to calculate the shortest path for all pairs of vertices in the graph with positive and negative weights. A naive method is to run the Bellman-Ford algorithm |V| times (|V| is the number of vertices). However, we'll see that another algorithm, combining SP and Bellman-Ford, can accomplish this task with better complexity.

## Problem definition

We have constructed a mountain village with N houses connected by forest roads, with varying elevations. The elevations are represented positively or negatively depending on the direction, the cost of going up being at least as large as the absolute value of the "negative cost" of going down. This means that we recharge our vehicle batteries during descents, while energy consumption is required during ascents. Here's an example of going downhill from house A to house B and uphill from B to A.

Example: A --(-2)--> B and B --(+7)--> A

The City Hall has a budget to buy K bread ovens, which can be installed in a village house (the ovens can't be placed outside houses, as they need to be supplied with electricity). The question is: how to place the K ovens to minimize the cost for the villager furthest from an oven?

We therefore have negative weights and, in addition, we have to calculate all the paths between two houses to determine where to place the ovens.

The village is represented as a weighted graph with positive or negative costs. It should be noted that there is no absorbing circuits in the village, as any upward movement is more costly than the inverse of any downward movement.

To solve this placement problem, we'll first revise the SP algorithm, then implement the new algorithm combining the SP and Bellman-Ford algorithms to calculate all the shortest paths between any two houses.

# Analyzing the complexity of the SP algorithm

In Lecture 2, you saw the complexity of the SP algorithm, depending on the data structure for the boundary (list or binary heap). Below you'll find the implementation with two data structures (SP_naive and SP_heap).

Give the worst-case time complexity of these two functions for sparse graphs {$|E|\approx|V|$} and dense graphs {$|E|\approx|V|^2$}, respectively. To do this, recall the complexity seen in Lecture 2.

def SP_naive (graph, s): frontier = [s] dist = {s:0} while len(frontier) > 0: x = min(frontier, key = lambda k: dist[k]) frontier.remove(x) for y, dxy in graph[x].items(): dy = dist[x] + dxy if y not in dist: frontier.append(y) dist[y] = dy elif dist[y] > dy: dist[y] = dy return dist from heapq import * def SP_heap (graph, s): frontier = [] heappush(frontier, (0, s)) done = set() dist = {s: 0} while len(frontier) > 0: dx, x = heappop(frontier) if x in done: continue done.add(x) for y, dxy in graph[x].items(): dy = dx + dxy if y not in dist or dist[y] > dy: heappush(frontier,(dy, y)) dist[y] = dy return dist

You now have functions that randomly generate sparse {$|E|\approx|V|$} and dense {$|E|\approx|V|^2$} graphs, as well as benchmark.

import random def random_sparse_graph (n, step): graph = {f'{i}_{j}': {} for i in range(1, n+1) for j in range(1, n+1)} for i in range(1, n+1): for j in range(1, n): d = random.randint(step+1, 2*step) graph[f'{i}_{j}'][f'{i}_{j+1}'] = d graph[f'{i}_{j+1}'][f'{i}_{j}'] = d for i in range(1, n): for j in range(1, n+1): d = random.randint(step+1, 2*step) graph[f'{i}_{j}'][f'{i+1}_{j}'] = d graph[f'{i+1}_{j}'][f'{i}_{j}'] = d return graph def random_dense_graph (n, d_max): graph = {f'{i}':{} for i in range(n)} for n1 in graph: for n2 in graph: if n2!= n1 and n2 not in graph[n1]: d = random.randint(1, d_max) graph[n1][n2] = d graph[n2][n1] = d return graph import matplotlib.pyplot as plt from timeit import timeit def benchmark(): Time_heap_sparse = [] Time_naive_sparse = [] Time_heap_dense = [] Time_naive_dense = [] n_list = [] for N in range(10, 30): n = N*N n_list.append(n) # on compare sur des cartes non-denses: grille de côté N contient N*N villes graph_sparse = random_sparse_graph(n = N, step = 100) # on calcule une moyenne sur N lancements en tirant aléatoirement une ville de départ à chaque fois Time_naive_sparse.append(timeit(lambda: SP_naive(graph_sparse, random.choice(list(graph_sparse))), number=N) / N) Time_heap_sparse.append(timeit(lambda: SP_heap(graph_sparse, random.choice(list(graph_sparse))), number=N) / N) # on compare sur des cartes denses graph_dense = random_dense_graph(n = N*N, d_max = 10000) Time_naive_dense.append(timeit(lambda: SP_naive(graph_dense, random.choice(list(graph_dense))), number=N) / N) Time_heap_dense.append(timeit(lambda: SP_heap(graph_dense, random.choice(list(graph_dense))), number=N) / N) plt.xlabel('N') plt.ylabel('T') plt.plot(n_list, Time_naive_sparse, 'r^', label="naive sparse") plt.plot(n_list, Time_heap_sparse, 'b^', label="heap sparse") plt.plot(n_list, Time_naive_dense, 'r*', label="naive dense") plt.plot(n_list, Time_heap_dense, 'b*', label="heap dense") plt.legend() plt.show()

Run the benchemark and explain the results by answering the following questions:

- According to your answer to the previous question, which structure has the better complexity for dense graphs in theory? What do you observe in the benchmark result? Can you explain?
- For sparse graphs, does the running result of the benchmark for the naive structure (list) correspond to the theoretical time complexity? Can you explain this?
- What is your general conclusion?

## You cannot run the benchmark?

You can click ici to see the result.

# Rewriting technique for new non-negative weights while preserving shortest paths

Given a graph with a weight function {$w$} that may have a negative-weight edge but no negative cycle, the idea of the rewriting technique is to compute a new set of non-negative weights while preserving the shortest paths so that Dijkstra's algorithm can be directly applied. More precisely, we need to find a new weight function {$\hat{w}$} such that the following two properties are satisfied:

- shortest path preservation: if {$p=\langle v_0, v_1, ..., v_k \rangle$} is any path from {$v_0$} to {$v_k$}, then {$w(p)$} is the shortest path from {$v_0$} to {$v_k$} if {$\hat{w}(p)$} is the shortest path from {$v_0$} to {$v_k$}.
- non-negative weights: for all edges {$(u, v)$}, the new weight {$\hat{w}(u, v)$} is non-negative.

These two properties can be elegantly obtained using the Bellman-Ford algorithm. We'll proceed step by step in what follows.

Given a graph {$G=(V, E)$} with weight function {$w$}, we construct a new graph {$G'=(V', E')$} with weight function {$w'$} such that

-{$V'=V\cup\{s\}$}, and {$s\notin V$}

-{$E'=E\cup\{(s, v) : v\in V\}$}

-{$w'(v, v')=0$} when {$v=s$}, {$w'(v, v')=w(v, v')$} otherwise.

Complete the function {$add\_source(graph, src)$} which takes a graph {$graph$} and a source vertex {$src$} as parameters and generates a new graph as described above.

def add_source(graph, src): if src in graph: return "error: the source vertex is already in the graph" graph_src=graph.copy() ############TODO : complete code######### return graph_src

Test your function with the following code:

sim_graph = {"A":{"B":-5,"C":2,"D":3}, "B":{"A":6,"C":4}, "C":{"D":1}, "D":{}} src="source" sim_graph_src=add_source(sim_graph, src) assert sim_graph_src == {"A":{"B":-5,"C":2,"D":3}, "B":{"A":6,"C":4}, "C":{"D":1}, "D":{}, "source":{"A":0, "B":0, "C":0, "D":0}}

Now complete the following function {$bellman\_ford(graph, src)$} which takes {$graph$} and {$src$} as parameters and returns a dictionary when there are no negative cycle: a key is a vertex in {$graph$} whose value is the distance of the shortest path from {$src$} to that vertex. It returns "None" when there is a negative cycle. You saw the Bellman-Ford algorithm in the course, which sets a target vertex. You'll use the same principle when setting a source vertex: the formula becomes {$\mbox{OPT}(i,v) = \min\Big(\mbox{OPT}(i-1,v), \min_{u\in V} (\mbox{OPT}(i-1,u) + \omega((u,v)))\Big)$}, where {$OPT(i,v)$} is the length of the shortest path from the source vertex to a vertex {$v$} that contains at most {$i$} arcs.

import math def bellman_ford(graph, src): dist={} ############TODO : complete code############# return dist

Test your function with the following code.

opt_distance=bellman_ford(sim_graph_src, src) assert opt_distance=={'A': 0, 'B': -5, 'C': -1, 'D': 0, 'source': 0} neg_cycle_graph = {"A":{"B":-5,"C":2,"D":3}, "B":{"A":3,"C":4}, "C":{"D":1}, "D":{}} assert bellman_ford(neg_cycle_graph, "A")==None

Complete the {$rewrite\_weights(graph, dist)$} function with two parameters: {$graph$} and {$dist$}, {$dist$} being a dictionary like the result returned by the {$bellman\_ford$} function. What the {$rewrite\_weights(graph, dist)$} function returns is a new graph calculated from the {$graph$} parameter as follows: for each {$graph[i][j]$} weight, the new value is obtained by {$graph[i][j]+dist[i]-dist[j]$}.

import copy def rewrite_weights(graph, dist): # use deepcopy altered_graph = copy.deepcopy(graph) # Recalculate the new nonnegative weights ############### TODO : complete code ################## return altered_graph

Test your function with the following code.

opt_distance={'A': 0, 'B': -5, 'C': -1, 'D': 0, 'source': 0} nonneg_graph=rewrite_weights(sim_graph, opt_distance) assert nonneg_graph=={'A': {'B': 0, 'C': 3, 'D': 3}, 'B': {'A': 1, 'C': 0}, 'C': {'D': 0}, 'D': {}}

Bonus: Why is the new weight ({$graph[i][j]+dist[i]-dist[j]$}) always non-negative? Can you demonstrate this?

# Calculation of shortest paths for all pairs

Complete the {$all\_distances(graph)$} function to calculate the shortest distance of each pair of vertices in {$graph$} (this graph contains only non-negative weights, e.g. the one obtained in the previous question) by calling the {$SP\_heap$} function for each vertex. What the {$all\_distances(graph)$} function returns is a dictionary whose key is a tuple of two vertices {$(v1, v2)$} and whose value is the shortest distance from v1 to v2. The value is None if v2 cannot be reached from v1.

def all_distances(graph): d = {(u,v):None for u in graph for v in graph} ############### TODO : complete code ################## return d

Test your function with the following code.

nonneg_graph={'A': {'B': 0, 'C': 3, 'D': 3}, 'B': {'A': 1, 'C': 0}, 'C': {'D': 0}, 'D': {}} assert all_distances(nonneg_graph)=={ ('A', 'A'): 0, ('A', 'B'): 0, ('A', 'C'): 0, ('A', 'D'): 0, ('B', 'A'): 1, ('B', 'B'): 0, ('B', 'C'): 0, ('B', 'D'): 0, ('C', 'A'): None, ('C', 'B'): None, ('C', 'C'): 0, ('C', 'D'): 0, ('D', 'A'): None, ('D', 'B'): None, ('D', 'C'): None, ('D', 'D'): 0 }

From all the above functions, now implement a {$BF\_SP\_all\_pairs(graph)$} function to calculate the shortest paths of all pairs for a graph that can contain negative weights. Please note that the expected results are the shortest paths for the original graph (with negative weights), not those for the graph with only non-negative weights obtained after rewriting.

Test your function with the following code.

assert BF_SP_all_pairs(sim_graph)== {('A', 'A'): 0, ('A', 'B'): -5, ('A', 'C'): -1, ('A', 'D'): 0, ('B', 'A'): 6, ('B', 'B'): 0, ('B', 'C'): 4, ('B', 'D'): 5, ('C', 'A'): None, ('C', 'B'): None, ('C', 'C'): 0, ('C', 'D'): 1, ('D', 'A'): None, ('D', 'B'): None, ('D', 'C'): None, ('D', 'D'): 0}

## aide

- The source vertex should be extracted from the graph before rewriting the weight of the graph.
- What you calculate with the function {$all\_distances$} are all shortest paths for the graph after rewriting. You need then to undo the rewriting to get the all shortest paths for the original graph: you can recover the original distances from the weights obtained by the Bellman-Ford algorithm.

What is the complexity of your algorithm in the previous question? Is there a difference when the input graph is sparse or dense? Please explain.

# Application (bonus): Placement problem

Let's return to our placement problem, where we need to calculate the shortest paths for all pairs of a given set of houses. Specifically, we need to write a program that allows the City Hall to choose the houses in which bread ovens can be installed, in order to minimize the cost of access to an oven for the villagers.

The input data are :

- a set of houses connected by forest roads, whose weight can be positive or negative;
- a subset of candidate houses to host a bread oven (it may be the case that not all houses can host an oven).

The program generates a subset of candidate houses in such a way that every house in the village has an oven or access to another house's oven. The idea is to take into account the cost of accessing the nearest house with an oven.

## Data representation

To represent the map of houses in the village, we'll use a Python **dictionary**.

Each house will be represented by a dictionary entry whose key is the name of the house (a string) and whose value is another dictionary whose keys are the names of neighboring houses and whose values are the access costs to these neighbors:

Consequently, a map can be presented as follows:

village = { 'Birchwood':{'Corner':15, 'Fairview':-6, 'Gables':8}, 'Corner':{...}, ... }

We see two nested dictionaries. The use of dictionaries allows access time in $O(1)$.

The village information will be provided in a file consisting of two parts separated by a line containing two dashes: --. The first part contains the list of houses. The second part contains the list of roads between the houses and the access costs between them, for example :

Birchwood Corner Fairview Gables -- Birchwood:Corner:15 Birchwood:Fairview:-6 Corner:Fairview:8

To load a card from a file, you can use the following code:

def read_map(filename): f = open(file=filename, mode='r', encoding='utf-8') map = {} while True: # reading list of cities from file ligne = f.readline().rstrip() if (ligne == '--'): break info = ligne.split(':') map[info[0]] = {} while True: # reading list of distances from file ligne = f.readline().rstrip() if (ligne == ''): break info = ligne.split(':') map[info[0]][info[1]] = int(info[2]) return map

Copy the toy example provided below (the toy_village dictionary) into your python file:

toy_village = {'A': {'B': -3, 'E': 20, 'F': 30}, 'B': {'A': 6, 'C': 9, 'F': 39}, 'C': {'B': 9, 'D': 8}, 'D': {'C': 8, 'F': 50}, 'E': {'A': -10, 'F': 6}, 'F': {'A': -20, 'B': -25, 'D': -15, 'E': 6} }

# Problem precision

Recall that the City Hall has funding to equip {$k$} houses for the bread oven. The number {$k$} is given in the problem instance.

For each house, we calculate the cost of access to the nearest equipped house (see Question 10). For a given set of houses, the access cost is the maximum cost of this set (see Question 11).

The aim is, for a village and a number {$k$} of bread ovens, to choose the location of {$k$} houses among the candidate houses in such a way as to minimize the overall access cost. Note that the number of candidate houses must not be less than {$k$}.

Write a function {$closest\_oven(house, oven\_houses, distance\_dict)$} taking as parameters a given house {$house$}, a list of houses equipped with an oven {$oven\_houses$}, and a dictionary of costs {$\{(u,v) : d(u,v) ...\}$} between all the houses, the latter being calculated in the previous section. The function returns a pair {$(d, v)$}, where {$d$} is the access cost (the distance between {$house$} and the nearest house in {$oven\_houses$}) and {$v$} is this nearest house equipped with an oven. If the list of houses with ovens given as an argument is empty, the function returns {$(math.inf, None)$}.

def closest_oven(house, oven_houses, distance_dict): ...

Test your function with the following code:

toy_distance_dict = BF_SP_all_pairs(toy_village) assert closest_oven('A', ['B','E','D'], toy_distance_dict) == (-3, 'B') assert closest_oven('B', ['B','E','D'], toy_distance_dict) == (0, 'B') assert closest_oven('C', ['B','E','D'], toy_distance_dict) == (8, 'D') assert closest_oven('D', ['B','E','D'], toy_distance_dict) == (0, 'D') assert closest_oven('E', ['B','E','D'], toy_distance_dict) == (-19, 'B') assert closest_oven('F', ['B','E','D'], toy_distance_dict) == (-25, 'B')

Write a function {$kcentre\_value(village, oven\_houses, distance\_dict)$}, which returns the maximum cost between a house in the village and the nearest oven, i.e. the access cost of this configuration of houses and houses equipped with ovens.

def kcentre_value(village, oven_houses, distance_dict): ...

Test your function using the following code:

assert kcentre_value(toy_village, ['B','E','D'], toy_distance_dict) == 8

To execute the following code, you must first download a village of 20 houses available here. Then understand the {$brute\_force$} function. What do you observe as the result of {$BF\_benchmark$}? Explain.

from itertools import * def brute_force(map, candidates, k, distance_dict) : best_combi = [] best_dist = math.inf for combi in combinations(candidates, k): dist= kcentre_value(map, list(combi), distance_dict) if dist<best_dist: best_combi= list(combi) best_dist= dist return best_dist, set(best_combi) assert brute_force(toy_village, list(toy_village), 3, toy_distance_dict) == (0, {'C', 'D', 'B'}) assert brute_force(toy_village, list(toy_village), 2, toy_distance_dict) == (8, {'C', 'A'}) import matplotlib.pyplot as plt from timeit import timeit village = read_map('village.map') village_distance_dict = BF_SP_all_pairs(village) def BF_benchmark(): Time_brute_force = [] k_list = [] for k in range(1, 20): Time_brute_force.append(timeit(lambda: brute_force(village, list(village), k, village_distance_dict), number=1)) k_list.append(k) #print N_list, time_list plt.xlabel('k') plt.ylabel('T') plt.plot(k_list, Time_brute_force, 'r^') plt.xticks(k_list) plt.show()

The itertools module provides a combinations function that allows you to generate all possible combinations.

# Going further: more other algorithms

We are now interested in a greedy algorithm (a dedicated heuristic) that would be both fast and intelligent enough to produce solutions that are not too far from the optimal solution.

To implement it, we build a list of {$k$} houses selected to host a bread oven. Initially, this list is empty. At each step, a new house is added to the list. The chosen house is the one which, in addition to the previously selected houses, minimizes the cost of access for the whole village.

Write a function {$greedy\_algorithm(map, candidates, k, distance\_dict)$} that builds a solution according to the greedy algorithm described above. The function returns a pair consisting of the access cost and the set of candidate houses chosen to host the ovens.

Test your function with the following code: compare the result with {$brute\_force$} on {$village.map$} for {$k=5$}?

force_d, force_h = brute_force(village, list(village), 5, village_distance_dict) greed_d, greed_h = greedy_algorithm(village, list(village), 5, village_distance_dict) assert greed_d >= force_d assert brute_force(village, list(village), 1, village_distance_dict) == (14, {'CL5'}) assert brute_force(village, list(village), 5, village_distance_dict) == (2, {'CL4', 'CL7', 'CL6', 'CL5', 'CL2'}) print(greed_d, force_d)

The next algorithm we are designing is a random selection algorithm for {$k$} houses among the candidate houses. This algorithm randomly tries a number of combinations of {$k$} candidate houses and retains the best solution among these combinations.

Write a function {$random\_algorithm$} that returns the best solution from a number of random solution trials:

def random_algorithm(map, candidates, k, distance_dict, trials=100) : ...

Here again, the function parameters are similar to those of the previous functions, with an additional parameter, trials, which designates the number of randomly generated solutions.

Now run the following benchmark. For all values of k on the village {$village.map$} in terms of execution time and solution quality. What are your conclusions?

def BF_G_R_benchmark(max_k, random_trials = 100): # compare on lite time_bruteforce_lite = [] time_random_lite = [] time_greedy_lite = [] d_bruteforce_lite = [] d_random_lite = [] d_greedy_lite = [] k_list = [] for k in range(2, max_k): k_list.append(k) time_bruteforce_lite.append(timeit(lambda: d_bruteforce_lite.append(brute_force(village, list(village), k, village_distance_dict)[0]), number=1)) time_greedy_lite.append(timeit(lambda: d_greedy_lite.append(greedy_algorithm(village, list(village), k, village_distance_dict)[0]), number=1)) time_random_lite.append(timeit(lambda: d_random_lite.append(random_algorithm(village, list(village), k, village_distance_dict, random_trials)[0]), number=1)) plt.subplot(2, 1, 1) plt.xlabel('k') plt.ylabel('T') plt.xticks(k_list) plt.plot(k_list, time_bruteforce_lite, 'r*') plt.plot(k_list, time_greedy_lite, 'b*') plt.plot(k_list, time_random_lite, 'y*') plt.subplot(2, 1, 2) plt.xlabel('k') plt.ylabel('d') plt.xticks(k_list) plt.plot(k_list, d_bruteforce_lite, 'r-') plt.plot(k_list, d_greedy_lite, 'b-') plt.plot(k_list, d_random_lite, 'y-') plt.show()