Table of contents

First Laboratory Session

The two laboratory sessions sessions of the Algorithmics and Complexity module are a small programming project. The problem to be solved is a land planning for ensuring access to a health care service. The goal is to place a certain number of hospitals in order to cover a territory at best.

Objectives of the TP

The work carried out during the two 3-hour practical sessions, completed possibly through personal work between the two sessions, has the following objectives:

1. To allow students to program classical algorithms seen in class:
   • graph search algorithm
   • shortest path algorithm
   • greedy algorithm
   • random algorithm
   • brute force algorithm (course 8)
2. To implement an environment and programming tools for the resolution of the problem:
   • reading entries
   • generation of random input sets
   • definition of metrics and realization of benchmarks.
   • display of inputs and outputs

Practical information

The development will be done under Python using an IDE with a debugger, such as Visual Studio Code that you used in the course of programming in the first sequence (SIP in SG1) and you learned to debug with it.

Definition of the problem

The objective of the project is to write a program that allows one to choose the city (or the cities) in which we set up hospitals to cover an entire area at best.

Input data is therefore :

• a set of cities linked by roads;
• a subset of candidate cities to host a hospital (not all cities can host a hospital).

The program has to produce a subset of candidate cities such as each city in the territory has either a hospital or is "covered" by a hospital of another city. This notion of coverage will be detailed in the second lab session. The idea is to take into account the distance between a city and the nearest hospital.

The program will also have a graphical user interface that allows us to display the result.

Data representation

We will define a data structure to represent our maps/cities in memory. This data structure must be adapted to be used in our algorithms and to be displayed in the graphical interface.

We will also need to import the map of the territory provided to the program as a file.

Map of the territory in memory

To represent the map of the territory in memory, we will use a Python dictionary.

Each city will be represented by a dictionary's entry whose key is the name of the city (a string) and the value is also a dictionary composed of two entries indexed by the following keys:

• coords : the Cartesian coordinates of the city on the map;
• next : the list of neighbors with their distances, which will again be a dictionary whose keys are the names of neighboring cities and the values are the distances (i.e. lengths of roads) to these cities.

Consequently, a map may be presented as:

map = {
       'La Baule-Escoublac':{'coords':(356, 638),
                            'next':{'Saint-Nazaire':15, 'Guérande':6, 'Saint-André-des-Eaux':8}},
       'Saint-Nazaire':{'coords':...},
       ...
}

Observe three nested dictionaries. The use of dictionaries (unlike lists) allows an optimal access time in {$\mathcal{O}(1)$}

Map of the territory as a file

The map of the territory will be provided to you in a file composed of two parts separated by a line containing two dashes: --. The first part contains the list of cities and their Cartesian coordinates. The second part of the file contains the list of direct connections between the cities and the distances of the roads between them, for instance:

La Baule-Escoublac:356:638
Saint-Nazaire:484:563
Guerande:180:340
Saint-André-des-Eaux:345:402
--
La Baule-Escoublac:Saint-Nazaire:15
La Baule-Escoublac:Guérande:6
La Baule-Escoublac:Saint-André-des-Eaux:8

To load a map 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]] = {'coords': (int(info[1]), int(info[2])), 'next':{}}

    while True:  # reading list of distances from file
        ligne = f.readline().rstrip()
        if (ligne == ''):
            break
        info = ligne.split(':')
        map[info[0]]['next'][info[1]] = int(info[2])
        map[info[1]]['next'][info[0]] = int(info[2])

    return map

Toy example

Create a tp.py file in your development environment and copy the toy example provided below (the toy_map dictionary) into it:

toy_map = {'A': {'coords': (100, 100), 'next': {'B': 140, 'E': 200, 'F': 200}}, 
           'B': {'coords': (200, 140), 'next': {'A': 140, 'C': 260, 'F': 140}}, 
           'C': {'coords': (360,  40), 'next': {'B': 260, 'D': 360}}, 
           'D': {'coords': (320, 360), 'next': {'C': 360, 'F': 280}}, 
           'E': {'coords': ( 80, 280), 'next': {'A': 200, 'F': 120}}, 
           'F': {'coords': (160, 240), 'next': {'A': 200, 'B': 140, 'D': 280, 'E': 120}} }

assert read_map('territoire_jouet.map') == toy_map

This territory corresponds to the data contained in the file territoire_jouet.map available here. Download this file and save it in your working directory. Check that the map loaded from territoire_jouet.map corresponds to toy_map.

Graphical visualization

We keep the tkinter library that you used during the programming course in SG1. To display a map in a graphical window, we provide you with the function draw_map(map, hospitals = []) which takes as parameters :

 - map : the dictionary representing the map in the memory format detailed in the previous section ;
 - hospitals: a Python list of the names of the cities hosting a hospital. These cities will be colored in red. This parameter is optional ;
 - l and h : the width and height of the display window, which should be sufficient to display all cities whose coordinates correspond to pixels.

Test the following code:

import tkinter as tk

def draw_map(map, hospitals = [], l = 800, h = 800):
    w = tk.Tk()
    c = tk.Canvas(w, width=l, height=h)

    for city in map:
        for vname in map[city]['next']:
            if city < vname:
                x1 , y1 = map[city]['coords']
                x2 , y2 = map[vname]['coords']
                c.create_line(x1, y1, x2, y2, dash=True)
                c.create_text((x1+x2)//2, (y1+y2)//2, text=map[city]['next'][vname])

    for city in map:
        x , y = map[city]['coords']
        color = 'red' if city in hospitals else 'black' 
        c.create_rectangle(x-5, y-5, x+5, y+5, fill=color)
        c.create_text(x+12, y+12, text=city)

    c.pack()
    w.mainloop()

draw_map(toy_map, ['A', 'C']) # remember to comment this line after the test

The result obtained should be as follows:

Automatic generation of maps

In order to be able to feed the implemented algorithms with data and make benchmarks to compare their performance, we need maps of different sizes (i.e. number of cities) and different densities (i.e. number of roads).

The function random_grid_map(n, step) takes as parameters two integers n and step and returns a square grid of {$n \times n$} cities. The coordinates of the cities (the unit is the pixel) are spaced by the step step% and each edge of the grid is labelled with a distance (an integer value) randomly drawn in the open interval {$]step, 2 \times step [$}.

import random

def random_grid_map (n, step):
    map = {f'{i}_{j}': {'coords':(j*step,i*step), 'next':{}} 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)
            map[f'{i}_{j}']['next'][f'{i}_{j+1}'] = d
            map[f'{i}_{j+1}']['next'][f'{i}_{j}'] = d

    for i in range(1, n):
        for j in range(1, n+1):
            d = random.randint(step+1, 2*step)
            map[f'{i}_{j}']['next'][f'{i+1}_{j}'] = d
            map[f'{i+1}_{j}']['next'][f'{i}_{j}'] = d

    return map

draw_map(random_grid_map(5,100))

For example, an execution of random_grid_map (5, 100) should generate a map as displayed below, only distances can vary :

You can use the following check_map(map) function to verify that your generated map is well formed:

def check_map(map, complete=False):
    for city in map:
        for neighbor in map [city]['next']:
            if neighbor not in map:
                print('Error: The neighbor', neighbor, 'of the city', city, 'does not exist')
                return False
            if city not in map [neighbor]['next']:
                print('Error: The city', city, 'and the city', neighbor, 'are not reciprocally neighbors')
                return False
            if map [city]['next'] [neighbor] != map [neighbor]['next'] [city]:
                print('Error: The city', city, 'and the city', neighbor, 'are not at the same distance reciprocally')
                return False
        if complete:
            if len(map[city]['next']) < len(map) - 1:
                print('Error : The graph is not complete, there are missing neighbors to :', city)
                return False
    return True

def Q_1():
    complete_map = random_dense_map(10,500)
    assert len(complete_map) == 10 and check_map(complete_map, complete=True) 
    draw_map(complete_map)

Computing the shortest paths

In the next part of the topic, we will need to compute the shortest distances between cities. We will use the algorithm of the shortest paths from a source s. This algorithm was invented by one of the fathers of computer science: Edsger Dijkstra.

We are only interested in calculating the distances between s and other cities: the path is not important for our problem.

Below is the Python version of the algorithm as presented in lecture 2, in which we have removed the parent dictionary. PCCs stands for Plus Courts Chemins in french.

def PCCs(graph,s):
    frontier = [s]
    dist = {s: 0}

    while len(frontier)>0:

        ### extraction of the node of the boundary having the minimum distance ###
        x = extract_min_dist(frontier,dist)

        ### for each updated neighbor of the border and its distance ###
        for y in neighbors(graph, x):
            if y not in dist:
                frontier.append(y)

            new_dist = dist[x] + distance(graph,x,y)
            if y not in dist or dist[y] > new_dist:
                dist[y] = new_dist

    return dist

In lecture 2 we discussed the impact of the data structure chosen to represent the frontier on the algorithm performance. We remind you that the complexity of such an algorithm is {$\mathcal{O}(|V|\times C_{extract\_min}\,+\,|E|\times C_{update})$}.

In the version above, this is the naive implementation where we use an unsorted Python list to represent the frontier. We want to compare it with an implementation where we use a binary heap ({$heap_{min}$}) to represent the frontier. Remember that each implementation has different {$C_{extract\_min}$} and {$C_{update}$} complexities. We also remind that the choice of data structure determines whether the worst case complexity of the algorithm depend or not on the density of the graph. The following table summarizes the set of worst-case complexities already analyzed in lecture 2.

We propose to verify experimentally these worst case complexities inputting randomly generated maps with different densities.

Naive implementation.

def Q_2():
    assert PCCs_naive(toy_map, 'A') == {'A': 0, 'B': 140, 'E': 200, 'F': 200, 'C': 400, 'D': 480}

Binary heap

The standard Python library includes the module heapq (documentation here) which provides us with the heap data structure. This library is installed by default.

After reading the documentation, play with the following example :

from heapq import *

tas = [(1,'a'), (2,'b')]      ### the french word for heap is 'tas'         
heapify(tas)               ### create a min heap from any list in O(n)
heappush(tas, (0,'c'))      ### inserting in O(log n) 
print(heappop(tas))        ### extraction of (0, 'c') in O(log n)
print(heappop(tas)[0])     ### extraction >>> '1'.
print(len(tas))            ### measure in O(1) >>> 1 
print(heappop(tas)[1])     ### extraction >>> 'b'

def Q_3():
    assert PCCs_heap(toy_map, 'A') == {'A': 0, 'B': 140, 'E': 200, 'F': 200, 'C': 400, 'D': 480}

Calculation of all the shortest paths

assert all_distances(toy_map) == {('A', 'A'):   0, ('A', 'B'): 140, ('A', 'C'): 400, ('A', 'D'): 480, ('A', 'E'): 200, ('A', 'F'): 200, 
                                  ('B', 'A'): 140, ('B', 'B'):   0, ('B', 'C'): 260, ('B', 'D'): 420, ('B', 'E'): 260, ('B', 'F'): 140, 
                                  ('C', 'A'): 400, ('C', 'B'): 260, ('C', 'C'):   0, ('C', 'D'): 360, ('C', 'E'): 520, ('C', 'F'): 400,
                                  ('D', 'A'): 480, ('D', 'B'): 420, ('D', 'C'): 360, ('D', 'D'):   0, ('D', 'E'): 400, ('D', 'F'): 280,
                                  ('E', 'A'): 200, ('E', 'B'): 260, ('E', 'C'): 520, ('E', 'D'): 400, ('E', 'E'):   0, ('E', 'F'): 120, 
                                  ('F', 'A'): 200, ('F', 'B'): 140, ('F', 'C'): 400, ('F', 'D'): 280, ('F', 'E'): 120, ('F', 'F'):   0}

map_grid = random_grid_map(10,100)
assert all_distances(map_grid, PCCs_heap) == all_distances(map_grid, PCCs_naive)

map_dense = random_dense_map(100,1000)
assert all_distances(map_dense, PCCs_heap) == all_distances(map_dense, PCCs_naive)

Benchmark

The Python module timeit allows one to measure the execution time of a Python function passed as a parameter and launched number times.

Using the following code, compare the execution time of the two versions {$naive$} and {$heap$} of the shortest paths computations.

import matplotlib.pyplot as plt
from timeit import timeit

def benchmark():
    Time_heap_grid = []
    Time_naive_grid = []
    Time_heap_dense = []
    Time_naive_dense = []

    n_list = []

    for N in range(10, 25):
        print(N) # to follow its running

        n = N*N
        n_list.append(n)

        # we compare on non-dense maps: grid of side N contains N*N cities
        map_grid = random_grid_map(n = N, step = 100)

        # an average is calculated over N launches by randomly selecting a starting city each time
        Time_naive_grid.append(timeit(lambda: PCCs_naive(map_grid, random.choice(list(map_grid))), number=N) / N)
        Time_heap_grid.append(timeit(lambda: PCCs_heap(map_grid, random.choice(list(map_grid))), number=N) / N)

        # we compare on dense maps
        map_dense = random_dense_map(n = N*N, d_max = 10000)

        Time_naive_dense.append(timeit(lambda: PCCs_naive(map_dense, random.choice(list(map_dense))), number=N) / N)
        Time_heap_dense.append(timeit(lambda: PCCs_heap(map_dense, random.choice(list(map_dense))), number=N) / N)

    plt.subplot(2, 1, 1)
    plt.xlabel('N')
    plt.ylabel('T')
    plt.plot(n_list, Time_naive_grid, 'r^', label='naive grid')
    plt.plot(n_list, Time_heap_grid, 'b^', label='heap grid')
    plt.legend()

    plt.subplot(2, 1, 2)
    plt.xlabel('N')
    plt.ylabel('T')
    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()

List of candidate cities

Do not forget that our goal is to write a program that allows you to choose the city or the cities where hospitals should be located to cover the territory. From a description of the territory and a list of candidate cities, the program must select a subset of the candidate cities such that each city in the territory is either hosting a hospital, or covered by a hospital of another city.

To do this, we will now write a function which lists all the possible solutions, i.e. all the possible solutions are combinations of candidate cities. This function will be used in the second lab session in the implementation of the backtracking and approximation algorithms.

A combination is a sub-list of candidate cities. For example, if our candidates are ['A', 'B', 'C', 'D'], we have {$2^4 = 16$} combinations to cover the territory:

[]
['A']
['B']
['A', 'B']
['C']
['A', 'C']
['B', 'C']
['A', 'B', 'C']
['D']
['A', 'D']
['B', 'D']
['A', 'B', 'D']
['C', 'D']
['A', 'C', 'D']
['B', 'C', 'D']
['A', 'B', 'C', 'D' ]

Generators

We want to iterate with a loop for on the combinations of candidate cities without going through a list of exponential size. We should transform the function all_combis(candidates) into a generator.

Without going into too much detail, this is possible with the use of the keyword yield which allows one to "resume" the execution of the function where it was left at each iteration of the loop for.

Replace the instruction print in all_combis(candidates) with a yield to obtain a generator that can be used in the following way :

for combi in all_combis(candidates):
    ...

Only {$k$} hospitals can be built with the available funding. For this reason we will generate the combinations of size {$k$} only.

for combi in all_combis_k(['A','B','C','D'], k=2):
    print(combi)

fake_candidates = [f'city{i}''' for i in range(50)]
for combi in all_combis_k(fake_candidates, k=1):
    print(combi)

To go further...

To overcome the problem you have just observed, we will change the algorithm for enumerating all the combinations and we will use an implicit graph. An implicit graph is a graph of which we build the necessary parts as the search continues. As we will see later, this technique is often used when exploring state spaces (see lecture 8).

Neighborhood of a combination

We are going to represent the space of the combinations of the candidate cities by an oriented graph in which each node is a pair of lists:

• a list of the cities already chosen and which corresponds to the beginning of a possible combination
• a list of cities that are still candidates

To create the neighborhood of each node, we extend the list of selected cities with an additional candidate city from the second list. A node thus has as many neighbors as there are cities in its second list of candidate cities.

For example, the node (['A','C'], ['B','D','E']) represents the combination of size two which chooses 'A' and 'C' and has three neighboring nodes where each adds one city among the still candidate cities ['B','D','E'].

for v in neighbors((['A', 'C'], ['B', 'D', 'E']):
    print(v)

>>> (['A', 'C', 'E'], ['B', 'D'])
     (['A', 'C', 'D'], ['B'])
     (['A', 'C', 'B'], []) 

Exploring the combination tree

We now want to write a generator combi_gen(candidates) which takes as an argument a list of candidate cities and displays all possible combinations of candidate cities. Assume that there are no duplicates in your list.

Efficient generation of combinations of size {$k$}

['A', 'B', 'C']
['A', 'B', 'D']
['A', 'C', 'D']
['B', 'C', 'D']