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).

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()

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.

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}}

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

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': {}}

Calculation of shortest paths for all pairs

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
 }

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}

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$}.

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')

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

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()

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.

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)

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

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()