Algorithmic Modelling
Workshop

Tutor guide

Version	Date	Designer	Proofreader	Comments
3.0	15/03/2022	Benjamin COHEN BOULAKIA	Jean-François DOLLINGER	Code optimization : Use of tuples for the adjacency list Addition of a new implementation of the adjacency list with Dictionary and Set Performance comparison of the two implementations of the adjacency list Added questions on implementation optimization Addition of theoretical and technical questions (in particular on the multigraph aspect of the 7 bridges of Königsberg) Complete redesign of the algorithm code Added a calculation progress bar for experimental complexity A few more detailed explanations
2.0	16/04/2021	Benjamin COHEN BOULAKIA		Addition of the calculation of experimental complexity Redesign and addition of theoretical questions
1.0	29/05/2019	Sébastien RAMON	Benjamin COHEN BOULAKIA	Original version.

Workshop skills

• Modelling the data of a specific problem in the form of a graph.
• Presenting the implementations of the different possible representations of a graph (sagittal representation, adjacency list, adjacency matrix)
• Creating a program for finding the ‘optimal’ solution to a specific problem modelled by a graph
• Selecting the most efficient data structures in terms of computation time
• Determining and Analysing the experimental computation time of an algorithm

Goals

This Workshop will allow students to answer, in part, the problem of the prosit 1 of the course by programming in Python 3. It aims to represent graphs via adjacency lists and adjacency matrices, then to program algorithms for showing that a graph is Eulerian and providing an Eulerian cycle if it exists. The second part should lead them to think about what determines the speed of an algorithm, and especially to understand how the data size and organization impact performance. They must also begin to understand that studying these performances, whether experimentally or theoretically, is not at all trivial.

Philosophy and recommendations

The purpose of this workshop is not simply to code, but also to address the variants of data structures representing graphs, and the impacts that these choices will have in terms of development and performance simplicity. It is therefore important not to ignore the questions posed in the text, and to lead to a debate on these questions.

In addition, there are many graph libraries, thanks to which the user no longer has to worry about data structures. If students wonder about the usefulness of implementing such structures themselves, they should be made aware of the need to understand the principles on which these libraries are based, in order to be able to make the best use of them. The beginning of part 2. Representation of a graph via a adjacency matrix proposes a discussion on this aspect (the comparison of the different representations in terms of performance and memory footprint).

Finally, many links are listed throughout the workshop, students should be encouraged to visit them systematically. While some address scientific culture topics that go beyond the scope of the course, others provide important indications on the solution to be brought to certain problems.

Introduction

The goal is to provide the technical team with an ‘itinerary’ allowing them to pass through all the streets in the area exactly once, starting from the crossroad of rue Smatti with rue du Général Cohen Boulakia (denoted as 'Departure' on the diagram) and returning there. It is also necessary to determine whether it will be feasible to determine a similar route for a larger map, and at what cost.

But before working on this map, let's take the 7 bridges of Königsberg problem. This problem consists in determining whether or not there is a walking route through the streets of Königsberg that allows, from any given starting point, to cross each bridge exactly once, and to return to its starting point. This diagram can be translated into a graph where the areas of the city are vertices and where the bridges to be crossed are the edges of the graph.

In this graph (which is actually a multigraph), passing through all the bridges once is equivalent to passing through each edge only once. Can we apply this method to the map that Agathe gave us ?

Yes, you have to build a graph where the intersections of the streets in the area are vertices numbered from 1 to 11 and where the streets to be passed through are the edges.

Attention, in the sagittal representation, the vertices are named from 1 to 11. However, some of the structures that we will implement throughout the Workshop are in general indexed from 0. The students will have to be very vigilant not to forget the -1 or the +1 which will sprinkle their code.

If some students solve the problem by numbering their vertices from 0 onwards, insist that they keep this numbering, as it will make them aware of their algorithms’ side effects.

1. Representation of a graph via an adjacency list

1.1 Implementation of the adjacency list

That's all well and good, but we can't give these graphs in input parameter of the program in this sagittal form ! We will need to implement a data structure. We will represent our graph in memory by an adjacency list using the lists in Python.

Let's see the principle of the adjacency list on an example:

We have two lists :

• Head which contains the list heads
• Succ which contains the lists of successors for each list head.

But be careful, in this example, the graph is directed. We consider an undirected graph. How to do this ?

We only have to represent an edge $a~— b$ as two arcs : the arc $a \rightarrow b$ and the arc $b \leftarrow a$.

This means that an undirected graph will have the same adjacency list as the corresponding symmetrical directed graph.

Let's implement area A of Agathe's diagram according to this structure (using tuples rather than lists, since we won't change their content) :

head_zone_A = \
 (1,3,5,7,11,15,17,19,23,25,27,29) #SOLUTION
succ_zone_A = \
 (2,4,1,5,4,8,1,3,5,9,2,4,6,8,5,11,8,10,5,7,9,3,4,8,7,11,6,10) #SOLUTION

Same thing now with the adjacency list allowing to represent the 7 bridges of Königsberg problem. Initialize this adjacency list so that 'Head7Bridges' represents the list of the heads of lists and 'Succ7Ponts' represents the list of the lists of successors.

The vertices will be numbered as follows :

1. North will be vertex 1
2. Center will be vertex 2
3. East will be vertex 3
4. South_ will be the vertex 4.

But as we have seen, this graph is a bit special. What is this particularity ? How to represent this graph by adjacency list ?

The graph is not a simple graph but a multigraph. For example, there are two edges between North and West.

The representation by adjacency list works perfectly in this case. If there are two arcs $a \rightarrow b$, the vertex $b$ will appear twice in the list of successors of $a$.

Let's go ahead to represent this graph :

head_7_ponts = \
 (1,4,9,12,15) #SOLUTION
succ_7_ponts = \
 (2,2,3,1,1,3,4,4,1,2,4,2,2,3) #SOLUTION

1.2 Implementing the functions

In order to check that the data entered in these lists is correct, and to get a bit familiar with Python, we will implement two methods that take a graph as a parameter (in the form of the two lists representing the adjacency list) and a vertex (in the form of its number) :

• the degreSommetGraphe(matrix, head, succ, vertex) method, which returns the degree of the vertex
• the voisinsSommetGraphe(matrix, head, succ, vertex) method, which returns the list of neighbours of the vertex

And we're going to test that on the graph of area A.

In the function degreSommetGraphe(matrix,  head, succ, vertex), the parameter succ is actually useless. It was added so as not to give students a hint about the problem-solving.

Moreover, the case of a graph consisting of a single vertex without edge is particular. On the one hand, it is necessary to decide how to represent the list, on the other hand, to manage this specific case with regard to functions. Here, it is not very interesting, we assume that we will not encounter this situation.

def degreSommetGraphe(head, succ, sommet):
 degre =\
 head[sommet] - head[sommet-1] #SOLUTION
 return degree

def voisinsSommetGraphe(head, succ, sommet):
 liste_voisins =\
 succ[head[sommet-1]-1:head[sommet]-1] #SOLUTION
 return liste_voisins

print("### Degree of the vertices in the graph of the 7 Könisgberg bridges ###")
for sommet in range(1, len(head_7_ponts)):
 print("vertex", str(sommet), ":", str(degreSommetGraphe(head_7_ponts, succ_7_ponts, sommet)))
    
print("\n### Neighbours of the vertices in the graph of the 7 bridges of Könisgberg ###")
for sommet in range(1, len(head_7_ponts)):
 print("sommet", str(sommet), ":", str(voisinsSommetGraphe(head_7_ponts, succ_7_ponts, sommet)))

print("\n### Degree of the Area A graph vertices ###")
for sommet in range(1, len(head_zone_A)):
 print("vertex", str(sommet), ":", str(degreSommetGraphe(head_zone_A, succ_zone_A, sommet)))

print("\n### Neighbours of the Area A graph vertices ###")
for sommet in range(1, len(head_zone_A)):
 print("vertex", str(sommet), ":", str(voisinsSommetGraphe(head_zone_A, succ_zone_A, sommet)))

### Degré des sommets du graphe des 7 ponts de Könisgberg ###
sommet 1 : 3
sommet 2 : 5
sommet 3 : 3
sommet 4 : 3

### Voisins des sommets du graphe des 7 ponts de Könisgberg ###
sommet 1 : (2, 2, 3)
sommet 2 : (1, 1, 3, 4, 4)
sommet 3 : (1, 2, 4)
sommet 4 : (2, 2, 3)

### Degré des sommets du graphe de la Zone A ###
sommet 1 : 2
sommet 2 : 2
sommet 3 : 2
sommet 4 : 4
sommet 5 : 4
sommet 6 : 2
sommet 7 : 2
sommet 8 : 4
sommet 9 : 2
sommet 10 : 2
sommet 11 : 2

### Voisins des sommets du graphe de la Zone A ###
sommet 1 : (2, 4)
sommet 2 : (1, 5)
sommet 3 : (4, 8)
sommet 4 : (1, 3, 5, 9)
sommet 5 : (2, 4, 6, 8)
sommet 6 : (5, 11)
sommet 7 : (8, 10)
sommet 8 : (5, 7, 9, 3)
sommet 9 : (4, 8)
sommet 10 : (7, 11)
sommet 11 : (6, 10)

This representation by adjacency list seems very efficient to determine all the neighbours of a vertex. But is it as effective for all processing tasks as one could imagine on a graph ? What if the graph is directed ?

In fact, the access to the list of the neighbours of a vertex is fast, since the time is constant whatever the number of neighbours. A treatment on all the neighbours of a vertex will therefore be proportional to the number of neighbours, without any other factor coming into play, which is relatively efficient.

Checking the existence of an edge will also require going through all the neighbours of a vertex, but here it is not ideal, we can do better.

Another complex operation is adding or removing an edge. In fact, it will be necessary to perform two memory lags regarding a part of the list of neighbours (one for each end of the edge).

Finally, if we consider a directed graph, looking for the predecessors of a vertex is also a rather cumbersome operation.

Let's see if we can correct these defects.

1.3 Pythonesque Implementation

The implementation of the adjacency list that we have produced is very generic, it can be applied to most languages, with a minimum of adaptation. It is important to know how to implement. But we might also want to exploit Python’s characteristics a little more, especially its specialized data structures, in order to avoid the performance degradations at runtime, which we talked about. Specialists like to talk about a Pythonesque implementation (which follows the philosophy of Python, which exploits its possibilities to the fullest).

There are two aspects on which we could improve our implementation of adjacency list. Which ones ?

• We can improve the way a list of neighbours is represented
• We can improve the way a vertex is associated with this list

What do you propose as a solution for these two points ?

To associate a vertex with a datum, the most efficient way is to use a list of hashing. In Python, you can directly use the Dictionary, which implements a hashing function to associate a value with a key. Here, the value is the name of the vertex (its number in this case), and the value is the set of neighbours of the vertex.

For the set of vertices, depending on the usage, different structures can be used :

• If we often need to check the presence of a vertex in the neighbours of another vertex, the Set function is more appropriate, because it checks the existence of a value in $O(1)$ (attention, this structure is unsuitable for multigraphs as the instance of the 7 bridges of Könisgberg, because it doesn’t allow to store twice the same value).
• If we instead need to iterate through all the neighbours of a vertex :
  • If the graph changes over time, we can use a List
  • Otherwise, a Tuple
  In both cases, iterate on the set of neighbours will be done faster than with a Set

Be careful though, in this version, if you number the vertices from 1, you have to adapt the upper bound of the academic career :

for sommet in range(1, len(liste_zone_A)+1):

Implement the version you want, and reimplement the function to display the degree of the vertices. Keep the initial version, so we can compare their performance !

liste_zone_A =\
{ #SOLUTION in this version, we use Set
 1 : {2, 4}, #SOLUTION
 2 : {1, 5}, #SOLUTION
 3 : {4, 8}, #SOLUTION
 4 : {1, 3, 5, 9}, #SOLUTION
 5 : {2, 4, 6, 8}, #SOLUTION
 6 : {5, 11}, #SOLUTION
 7 : {8, 10}, #SOLUTION
 8 : {3, 5, 7, 9}, #SOLUTION
 9 : {4, 8}, #SOLUTION
 10: {7, 11}, #SOLUTION
 11: {6, 10} #SOLUTION
} #SOLUTION


def degreSommetGraphe2(liste, sommet):
 degre =\
 len(liste[sommet]) #SOLUTION
 return degree

print("### Degree of the Area A graph vertices ###")
for sommet in range(1, len(liste_zone_A)+1):
 print("vertex", str(sommet), ":", str(degreSommetGraphe2(liste_zone_A, sommet)))

### Degré des sommets du graphe de la Zone A ###
sommet 1 : 2
sommet 2 : 2
sommet 3 : 2
sommet 4 : 4
sommet 5 : 4
sommet 6 : 2
sommet 7 : 2
sommet 8 : 4
sommet 9 : 2
sommet 10 : 2
sommet 11 : 2

For the list of neighbours of the vertices, it’s even more direct because there’s no need to write a function, a simple print(graph) is enough !

But does it work with the 7 bridges of Königsberg problem ?

As with the generic version of the adjacency list, it works in the state. It’s not even possible to make mistakes when implementing functions.

And in terms of performance ? We will measure the computation time, with the magic command %%timeit. And we recall once again that :

$10^n$	Prefix	Symbol	Decimal
$10^{−1}$	deci	d	0.1
$10^{−2}$	centi	c	0.01
$10^{−3}$	milli	m	0.001
$10^{−6}$	micro	µ	0.000 001
$10^{−9}$	nano	n	0.000 000 001
$10^{−12}$	pico	p	0.000 000 000 001
$10^{−15}$	femto	f	0.000 000 000 000 001
$10^{−18}$	atto	a	0.000 000 000 000 000 001
$10^{−21}$	zepto	z	0.000 000 000 000 000 000 001
$10^{−24}$	yocto	y	0.000 000 000 000 000 000 000 001

%%timeit
# simple implementation
for sommet in range(1, len(head_zone_A)):
 res =\
 degreSommetGraphe(head_zone_A, succ_zone_A, sommet)

1.5 µs ± 11.5 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

%%timeit
# optimized implementation
for sommet in range(1, len(liste_zone_A)+1):
 res =\
 degreSommetGraphe2(liste_zone_A, sommet)

1.26 µs ± 21.6 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

It is indeed better in the optimized case, but we remain in the same orders of magnitude. That said, calculating the degree is one thing. But if we want to retrieve the list of the successors of a particular vertex ? Let's compare again :

%%timeit
# simple implementation
for sommet in range(1, len(head_zone_A)):
 indiceCourant =\
 head_zone_A[sommet-1] #SOLUTION
 indiceSuivant =\
 head_zone_A[sommet] #SOLUTION
 res =\
 succ_zone_A[indiceCourant-1:indiceSuivant-1] #SOLUTION

1.74 µs ± 33 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

%%timeit
# optimized implementation
for sommet in range(1, len(liste_zone_A)+1):
 res =\
 liste_zone_A[sommet] #SOLUTION

510 ns ± 61.4 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

This time, the Pythonesque version is much faster, we completely change the order of magnitude !

But there are operations that the adjacency list does not perform well, regardless of the implementation. Which ones ?

Finding the predecessors of a vertex in a directed graph remains a cumbersome operation. That said, we can compensate by also implementing an inverted adjacency list (a list of predecessors). This structure is sometimes called as adjacency map.

The conclusion is that the choice of implementation depends on the specific need you have. It's up to you to know the most typical and recurring use (systematic access, search, modification, etc.), in order to be able to select the most efficient structure.

2. Representation of a graph via a adjacency matrix

Since this representation by adjacency list has limits, let's try to use another structure : an adjacency matrix ! If we take the example from earlier, it gives this representation :

$$\begin{pmatrix} 0 & 1 & 0 & 0 & 0 & 0\\ 1 & 0 & 0 & 1 & 0 & 0\\ 1 & 1 & 0 & 0 & 1 & 0\\ 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1\\ 0 & 0 & 1 & 0 & 0 & 1 \end{pmatrix}$$

In this representation, we use a square matrix $m$ of size $n\times n$ (with $n$ the number of vertices) in which an edge $a~— b$ is represented by the value $1$ in the element $m_{ab}$. For vertices not connected to each other, we put the value $0$.

What do you think of this representation ? Is it better than the adjacency list ?

It is neither better nor worse. Some operations are easier to perform on this structure, and conversely some are easier to perform on the adjacency list
.

More specifically, the matrix is more efficient for accessing the predecessors of a neighbour (it’s as fast as accessing its successors), checking the existence of an edge, and modifying the graph (regarding the last two points, this difference is less prominent when using the Pythonesque version of the adjacency list).

On the other hand, building a list of neighbours of a vertex is more cumbersome, since it is necessary to browse the entire related line of the matrix, which implies a processing in $n$ operations ($n$ being the number of vertices).

There is another very important aspect for scientific computing, and which we have not mentioned so far - memory space. When we consider vertices that are counted in millions or even billions (which is the case of the social graph of Facebook, for example), it is essential. From this point of view, what is the most effective structure ?

We are tempted to answer that it is the adjacency list, since the matrix occupies in memory a place proportional to the square of the number of vertices squared, where the list occupies a place proportional to the number of edges. But in fact, these two values are close to each other when the graph is very dense. Since most operations on graphs are fast on one adjacency matrix (except for building a list of neighbours, but this is still acceptable), we tend to use the adjacency list for large low-density graphs, and the adjacency matrix for the rest.

Moreover, there are methods to represent a adjacency matrix more efficiently than the one seen here, in particular by using a hashing table indexed by pairs of vertices. This approach, although more memory intensive than the adjacency list, still has a smaller memory footprint than a classic matrix representation, and also favours the location of data in memory, which allows increased access performance in certain cases thanks to the memory paging mechanism.

2.1 Implementation of the adjacency matrix of Area A

We will start by initializing the adjacency matrix that allow us to represent area A of Agathe’s diagram : matrix_zone_A represents this matrix. Here, we do not consider distances, that is, the lengths of each street. Why ? What does this imply at the level of the matrix ? What structure should we use if we wanted to take these lengths into account ?

We might think that the sizes of the edges are important since in the long term, one wants to carry out a route which is as short as possible. But here it's actually useless, since we want to go through each edge only once. Whether or not the length of the edges is taken into account does not change anything. On the other hand, this will not be the case for the project
.

This is also another advantage of the adjacency matrix : representing values on the edges (this is known as a valued graph) is simpler than with an adjacency list, where the integers of the list of successors must be replaced by a two-integer structure (one for the number of the successor, one for the value of the corresponding edge, which in Python can be done with a tuple, or even better, a Named Tuple, or a dictionary if distances are subject to change).

Moreover, we consider an undirected graph. As a result, what structural property does this matrix have ?

It’s a symmetric matrix, since the existence of an edge $a~—~b$ will in fact be represented in our matrix by the arc $a~\rightarrow~b$ and the arc $b~\leftarrow~a$ (as with the adjacency list).

Let’s state this matrix. But this time, we are not going to use nested tuples, because the algorithm for calculating the Eulerian cycle that we are going to implement needs to modify the graph (without anticipating, it will gradually remove the edges already traversed from the graph). So, it will be a nested list.

matrix_zone_A = [
 [0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0], #SOLUTION
 [1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], #SOLUTION
 [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0], #SOLUTION
 [1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0], #SOLUTION
 [0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0], #SOLUTION
 [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1], #SOLUTION
 [0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0], #SOLUTION
 [0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0], #SOLUTION
 [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0], #SOLUTION
 [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1], #SOLUTION
 [0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0], #SOLUTION
]

In order to check that the data entered in this matrix is correct, we will implement the new versions of the previous functions, degreSommetGrapheMatrice(matrix, vertex) and voisinsSommetGrapheMatrice(matrix, vertex).

For the list of neighbours, the use of Comprehension lists will certainly be more practical than if we implemented a double loop nested (and faster).

For the degree, there are two different approaches. The first is to count the number of neighbours (with the previous function). There is another, a little more subtle. Do you see which one ?

We can also sum the values of the line corresponding to the vertex whose degree we want to calculate using the function sum()

There are two reasons to prefer this solution :

• The result is much faster than with the counting of neighbours by Comprehensions
• Unlike the Comprehension approach, this method also works with multigraphs, so with the graph of the 7 Königsberg bridges

First, let them choose the method they prefer. Even if they don't think about the case of multigraphs, they should have the reflex to test the performance of both implementations.

The choice is yours !

And be careful, as this time we use lists, the numbering of the vertices necessarily starts at 0. Let's still keep a consistent display with the version by adjacency list (it's a bad idea, nothing better to get tangled up, but in this case, we can compare the different versions).

Be careful with these advanced functions, as they can give the impression that building a list of neighbours is not expensive. Students should not fall into the trap. These approaches do not remove complexity from processing, they just make them more transparent when writing code. But the advantages and disadvantages of the data structures we manipulate are still there (even if these Python functions are optimized from the point of view of access performance).

Concerning the numbering displayed which differs from the internal indexing, it is indeed a bad idea, it is better to have a display consistent with the code. This situation illustrates well the fact that it is better to think before implement an algorithm or a data structure. And if the students decided from the start to number their vertices starting from 0, they were right !

def voisinsSommetGrapheMatrice(matrice, sommet):
 liste = matrice[sommet]
 voisins = \
 [i for i, value in enumerate(liste) if liste[i] != 0] #SOLUTION
 return voisins

def degreSommetGrapheMatrice(matrice, sommet):
 degre =\
 sum(matrice[sommet]) #SOLUTION
 return degree

print(“### Degree of vertices of the graph for Area A ###")
for sommet in range(len(matrix_zone_A)):
 print("vertex", str(sommet+1), ":", str(degreSommetGrapheMatrice(matrix_zone_A, sommet)))

print("\n### Neighbours of the Area A graph vertices ###")
for sommet in range(len(matrix_zone_A)):
 voisins = voisinsSommetGrapheMatrice(matrix_zone_A, sommet) # we proceed in two steps, because
 print("vertex", str(sommet+1), ":", str([v+1 for v in voisins])) # les indices commencent à 0

### Degré des sommets du graphe de la Zone A ###
sommet 1 : 2
sommet 2 : 2
sommet 3 : 2
sommet 4 : 4
sommet 5 : 4
sommet 6 : 2
sommet 7 : 2
sommet 8 : 4
sommet 9 : 2
sommet 10 : 2
sommet 11 : 2

### Voisins des sommets du graphe de la Zone A ###
sommet 1 : [2, 4]
sommet 2 : [1, 5]
sommet 3 : [4, 8]
sommet 4 : [1, 3, 5, 9]
sommet 5 : [2, 4, 6, 8]
sommet 6 : [5, 11]
sommet 7 : [8, 10]
sommet 8 : [3, 5, 7, 9]
sommet 9 : [4, 8]
sommet 10 : [7, 11]
sommet 11 : [6, 10]

2.2 Implementation of the adjacency matrix of the 7 bridges of Könisgberg problem

Same thing now with the adjacency matrix allowing to represent the 7 bridges of Königsberg problem. Initialize this matrix : matrix_7_ponts represents this matrix.

The vertices will be numbered as follows :

• North will be vertex 1
• Center will be vertex 2
• East will be vertex 3
• South will be vertex 4

As for the adjacency list, we must take into account the particularity of this graph. How can we represent this graph with a adjacency matrix ?

As we have already seen, we are dealing with a multigraph. In the representation by adjacency matrix, the solution consists in placing in the matrix the values representing the number of times the vertices are linked together. For example, the two edges connecting North and West will be represented by the value 2 in the matrix.

Fortunately, we don't need to store the lengths… But doesn't that pose a problem with respect to our functions ?

For the list of neighbours, you just have to make sure that the conditional Comprehension test is indeed if list[i] != 0 and not if list [i] == 1

For the degree of a vertex, we will not be able to count the number of non-zero values of the line corresponding to the vertex whose degree we want to calculate, we must use the version with the function sum(). Students who chose the wrong solution just have to recode this function ! Anyway, since the sum implementation was faster, they should have done the test and chosen this version...

matrix_7_ponts = [
 [0,2,1,0], #SOLUTION
 [2,0,1,2], #SOLUTION
 [1,1,0,1], #SOLUTION
 [0,2,1,0] #SOLUTION
]

print("### Degree of the vertices in the graph of the 7 Könisgberg bridges ###")
for sommet in range(len(matrix_7_ponts)):
 print(degreSommetGrapheMatrice(matrix_7_ponts, sommet))

### Degré des sommets du graphe des 7 ponts de Könisgberg ###
3
5
3
3

3. Eulerian Cycle

For the implementation of the algorithms of this part, we will use the adjacency matrix to facilitate the development. Indeed, as the algorithm goes through the edges, it will remove them from the graph, so it must be able to simply modify it (besides, be careful, once the algorithm is executed, we can no longer use the graph). This being the case, one could wonder if with the adjacency list implemented by dictionaries and set, it wouldn't be faster to execute.

3.1 Existence of an Eulerian cycle in a graph

Thanks to Euler’s theorem, we know that there is an Eulerian cycle in a graph if and only if the graph is connected and a specific condition on it is satisfied. What is this property that a connected graph must satisfy to be Eulerian ?

There is no odd degree vertex in the graph.

Complete the function existeCycleEulerien(matrix) which takes a connected graph as a parameter (in the form of an adjacency matrix) and returns True if there is an Eulerian cycle in the graph. You may well need the modulo !

Let’s go !

def existeCycleEulerien(matrice):
 for sommet in range(len(matrice)):
 degre = degreSommetGrapheMatrice(matrice, sommet) #SOLUTION
 if((degre % 2) > 0): #SOLUTION
 return False; #SOLUTION
 return True; #SOLUTION
  
if (existeCycleEulerien(matrix_zone_A)):
 print ("The graph of Area A is Eulerian")
else:
 print ("The graph of Area A is not Eulerian")

Le graphe de la Zone A est eulérien

Let’s also test the function on the graph of the 7 bridges of Könisgberg !

if (existeCycleEulerien(matrix_7_ponts)):
 print ("The graph of the 7 bridges of Könisgberg is Eulerian")
else:
 print ("The graph of the 7 bridges of Könisgberg is not Eulerian")

Le graphe des 7 ponts de Könisgberg n'est pas eulérien

It seems that Agathe was right !

3.2 Calculation of a Eulerian cycle in a graph

Now that we are able to determine if a graph is Eulerian, it is time to find a way to calculate a Eulerian cycle when we have such a graph.

This is exactly what the following algorithm does, using the principle of backtracking. Here is the algorithm :

function CycleEulérien(graphe)

Create a cycle and a stack of empty vertices
Initialize the current vertex as the first vertex of the matrix

Repeat until the stack is empty and the current vertex has no more neighbours :

If the current vertex has at least one neighbour :

Add the current vertex to the stack
We delete the edge between the current vertex and this neighbour
The current vertex becomes this neighbour

Otherwise :

we add the current vertex to the cycle (principle of Backtracking)
we remove the 1st element of the stack that becomes the current vertex

By the way, why use a stack ? What algorithmic method could have been used instead ?

Recursive functions. Moreover, algorithmically, these two approaches are strictly equivalent, since the execution contexts of the functions are stacked during the different calls. But from an implementation point of view, a stacking only of the data useful to the algorithm will always be more efficient than a complete context switch context switching managed by the OS.

In any case, we already know which data structure Python to use for this stack. Moreover, we will use the same one to store the cycle. Logically, we just add vertices at the end.

Complete the cycleEulerien(matrix) method, which uses the principle of this algorithm. It takes a graph with an Eulerian cycle as a parameter (in the form of a matrix) and returns one of the graph’s Eulerian cycles.

from collections import deque
import copy

def cycleEulerien(matrice):
 # the matrix is passed by reference, so we make a copy of the matrix to avoid overwriting its data.
 # as we must also copy the internal lists, we make a _deep copy_
 matrice = copy.deepcopy(matrice)
 n = len(matrice)

 cycle = deque() # cycle is the cycle to build
 stack = deque() # stack is the stack of vertices to process
 cur = 0 # cur is the current vertex. we start with the first vertex of the matrix

 # on loop as long as there are vertices to process in the stack
 # or the current vertex has at least 1 unprocessed neighbour
 while(len(stack) > 0 or degreSommetGrapheMatrice(matrice, cur) != 0):
          
 # if the current vertex has no neighbours, add it to the cycle
 # and return to the previously added vertex in the stack (backtracking)
 # which becomes the new current vertex
 if (degreSommetGrapheMatrice(matrice, cur) == 0): #SOLUTION
 cycle.append(cur) #SOLUTION
 cur = stack.pop() #SOLUTION
  
 # if it has at least 1 neighbour, we add it to the stack to come back to it later (backtracking)
 # we remove the edge it shares with this neighbour, which becomes the current vertex
 else:
 for i in range(n):
 if matrice[cur][i]: #SOLUTION
 stack.append(cur) #SOLUTION
 matrice[cur][i] = 0 #SOLUTION
 matrice[i][cur] = 0 #SOLUTION
 cur = i #SOLUTION
 break
 return cycle;

print("### Calculation of an Eulerian cycle of the graph of Area A ###")
cycle = cycleEulerien(matrix_zone_A)
for sommet in cycle:
 print(sommet+1, "-> ", end = '')
print(cycle[0]+1)

### Calcul d'un cycle Eulérien du graphe de la Zone A ###
1 -> 4 -> 9 -> 8 -> 7 -> 10 -> 11 -> 6 -> 5 -> 8 -> 3 -> 4 -> 5 -> 2 -> 1

Check that it is indeed an Eulerian cycle.

3.2. Computation time analysis

Our algorithm seems to work. And it is very fast. By the way, how long did it take to finish ? To check this, we could use the magic command %%timeit, or Python profilers like cProfile or profile. This kind of tool is very powerful to study the most time-consuming CPU or memory consuming code elements, their use is also strongly recommended if you need to optimize complex code.

But it would be interesting to display the evolution of the computation time according to the graph size. However, for that, it is necessary to store the times measured in variables to be able to then draw a curve, so we forget %%timeit and the profilers. The time library offers plenty of functions that might be useful. Which one do you think is suitable ?

The risk is that if we simply measure the time spent executing the function, so we will also include the pause times during which the OS does other things. The use of functions like time() is therefore excluded. The function process_time() is fine, because it only counts the CPU time spent on the current process.

If the students compare with the values displayed by %%timeit, they will see a small difference, nothing abnormal.

Let's test it right away !

import time

start = time.process_time() #SOLUTION
cycle = cycleEulerien(matrix_zone_A) #SOLUTION
stop = time.process_time() #SOLUTION

print(stop-start)

0.0

To be fast is fast ! But maybe on our instance, we were lucky, and another graph of the same size would be heavier to manage. How could we be a little more sure of the result ?

We would have to generate random graphs of the same size as the one we used, and calculate the average of computation times

We could do all this by hand, but there are libraries that will make our job easier. The NumPy library, a must-have for scientific computing in Python, could be useful to us. Among the Random sampling functions it offers, the function numpy.random.choice() allows, among other things, to generate a set of Booleans, distributed in the structure of our choice (here in a matrix), that we will only have to transform into integers.

Be careful all the same, we are considering an undirected graph, so we have to build a symmetric matrix. There is a simple technique : generating a Boolean matrix, and overlaying it with its transpose (a simple logical $or$ between the two matrices). All that’s left now is to convert the result into an integer matrix (in Python, the Boolean values True and False are represented by 0 and 1, respectively. If we had wanted to generate multigraphs, we would only have to add these two matrices instead of overlaying them. Rest assured, the code is provided to you !

There are still some problems with this approach. Do you see them ?

The resulting graph may not be Eulerian. Two approaches could remedy this problem :

• the rejection sampling. The general principle is complex, but here, it’s basically a matter of regenerating a new matrix until obtaining one that is Eulerian. This method is simple to implement, but it sometimes takes a long time to generate a valid matrix
• the design of an algorithm capable of generating a matrix with the right property. This approach remains very complicated (it’s at the level of a doctoral thesis) and is outside the scope of this Workshop

Note : students could begin to imagine a method in which we add edges between each pair of odd degree vertices if the graph is not connected. Nevertheless, this approach poses a problem : if the number of odd degree vertices is also odd, we will get an odd degree vertex at the end of these edge additions. Removing this vertex would create another odd degree vertex, as would removing an edge from this vertex. Especially since this method would pose the problem of the uniformity of the generation (see below).

Another problem is that the graphs we generate are not necessarily connected. This may also distort the results, since the cycle will only be on the connected subgraph in which the first vertex of the cycle is located, and the algorithm will end faster than expected. That said, modifying the algorithm so that it calculates an Eulerian cycle for each connected component of the graph is not very complicated (but it’s not the goal here).

Finally, there is no guarantee here that the graphs are generated in a uniform way, i.e. that they all have the same probability of being generated. The risk is to see certain graphs appear more often, and therefore to obtain biased results (this is why the methods proposed in the thesis mentioned above are so complex).

Given the complexity to solve all these problems, we will stick to the simple solution (without rejection or guarantee of connection) despite the defects it has. Nevertheless, this highlights the importance of being able to implement a high-performance uniform random generation algorithm if we wish to thoroughly study the experimental behaviour of an algorithm. It is clear that this problem requires the use of probabilistic notions that are sometimes quite complex (such as Markov chains).

You will find in the next cell the code generating this matrix. No need to understand it in detail (even if it is not very complicated thanks to NumPy), we will just use it.

import numpy as np

def grapheAleatoireEulerien(taille):
 # we generate a random Boolean matrix (to be able to do `or`)
 b = np.random.choice((True, False), size=(taille,taille), p=[0.4, 0.6])
    
 # we make the logical `or` of the matrix and its transpose
 b_symm = np.logical_or(b, b.T)
    
 # we return the Boolean matrix converted into integer matrix
 return b_symm.astype(int)

print(grapheAleatoireEulerien(4))

[[1 1 0 1]
 [1 0 1 1]
 [0 1 0 1]
 [1 1 1 0]]

Perfect ! Now we just have to test all this in loop on our algorithm. A hundred iterations should be enough.

import time

duree = 0
nb_iteration = 100
for i in range(nb_iteration):
 instance = grapheAleatoireEulerien(11) #SOLUTION
 start = time.process_time() #SOLUTION
 cycle = cycleEulerien(instance) #SOLUTION
 stop = time.process_time() #SOLUTION
 duree += (stop-start) #SOLUTION
print(duree/nb_iteration) #SOLUTION

0.00046875

Indeed, we see a slightly longer computation time appear, but nothing too worrying (reminder : the time is expressed in seconds). Besides, it seems logical !

But what would happen if the graphs were larger ? Let's try to vary this parameter, and see how the computation time evolves. We can display everything in the form of a curve. Matplotlib is another must-have for scientific computing, and it will be perfect for that. It’s a very powerful tool that can, for example, display a curve whose ordinate values are in a list. Let’s see this on an example :

from matplotlib import pyplot as plt
plt.plot([0, 0.301, 0.477, 0.602, 0.698, 0.778, 0.845, 0.903, 0.954, 1])
plt.show()

With this, we can display the evolution of the average computation time according to the size of graphs.

And then we can also determine the best-case and worst-case computation times. This will give us a slightly more detailed view of how the algorithm works. It’s all the more interesting to compare them because, although the worst-case complexity is still fairly easy to study (actually, this is one of the Prosit’s goals), the [average-case complexity]https://en.wikipedia.org/wiki/Average-case_complexity) is much more complex to address, and requires advanced mathematical tools, especially in enumeration.

And during the passing process, we will display a slightly larger curve, with the appropriate captions (be careful, if your computer lacks power, reduce the maximum size; you can also consider certain graph sizes only thanks to the step parameter of the range() function (nevertheless, keep values that are multiples of 10), but avoid reducing the number of iterations, otherwise you will have unrepresentative results).

Finally, as the calculations will start to take time, we will add a [progress bar] (https://ipywidgets.readthedocs.io/en/latest/examples/Widget%20List.html#IntProgress) to know where we stand. The code is provided to you (this time...).

from ipywidgets import IntProgress
from IPython.display import display

nb_iteration = 50
taille_min = 10
taille_max = 200
taille_step = 10

durees_moy = []
durees_mieux = []
durees_pire = []

# we display the progress bar
nb_tests = ((taille_max - taille_min) / taille_step) * nb_iteration
bar = IntProgress(min=0, max=nb_tests, layout={"width" : "100%"})
display(bar)

for taille in range(10, 200, 10):
 duree_mieux = float('inf') # yes, Python allows us to manipulate infinite values
 duree_pire = 0.0
 duree_moy = 0.0
    
 for i in range(nb_iteration):
 instance = grapheAleatoireEulerien(taille) #SOLUTION
 #SOLUTION
 start = time.process_time() #SOLUTION
 cycle = cycleEulerien(instance) #SOLUTION
 stop = time.process_time() #SOLUTION
 duree = (stop-start) #SOLUTION
 #SOLUTION
 duree_moy += duree #SOLUTION
 if duree_mieux > duree: #SOLUTION
 duree_mieux = duree #SOLUTION
 if duree_pire < duree: #SOLUTION
 duree_pire = duree #SOLUTION
            
 bar.value += 1 # we update the progress bar
        
 # we update the running time lists
 durees_moy.append(duree_moy/nb_iteration)
 durees_mieux.append(duree_mieux)
 durees_pire.append(duree_pire)

# we hide the progress bar
bar.close()

# we adjust the display of curves
tailles = [x for x in range(10, 200, 10)]
plt.figure(figsize=(20,10))

plt.xlabel('graph size')
plt.xticks(ticks=tailles) # valeurs affichées sur l'axe X
plt.ylabel('computation time')

# we load the data
plt.plot(tailles, durees_moy, label='average duration')
plt.plot(tailles, durees_mieux, label='best-case duration')
plt.plot(tailles, durees_pire, label='worst-case duration')

# we display
plt.legend()
plt.show()

IntProgress(value=0, layout=Layout(width='100%'), max=950)

What do you observe ? How can these results be interpreted ?

The results may take a few minutes to calculate. If some students have low capacity computers, display the correction graph.

In all three cases, we observe that the computation time increases faster than linearly. This means that multiplying the number of vertices by $n$ will do more than multiplying the computation time by $n$. The computation time grows faster than the number of vertices of the graph.

Moreover, the best-case and worst-case curves are much more 'uneven', which makes sense, since an average produces 'smoothed' values. If we increased the number of iterations and relied on a more uniform random generator, we might observe something more regular.

On the other hand, the three curves are quite close. This means that random graphs are all about equally difficult to process by our algorithm. In any case, and once again, on the graphs constructed by our random generator.

Nevertheless, the curve at best is sometimes closer to the curve on average than the curve at worst (for the size 160, for example). This is a case where most instances behave more or less the same, except for one that took a little longer.

It’s important to emphasize to students the need for keeping a large number of instances for each size, if a clear trend is to be obtained.

Could this result have been anticipated ?

Yes, by performing a theoretical study of the asymptotic temporal complexity of the algorithm :

There are $\lvert V\rvert$ edges to delete, and the list of vertices contains, at most, $\lvert E\rvert-1$ vertices in the list. The asymptotic complexity is therefore $O(\lvert V\rvert\times\lvert E\rvert)$. On average, it is much lower, because on the one hand each passage on a vertex decreases the number of edges connected to this vertex, and on the other hand $\lvert E\rvert-1$ is the maximum degree of a vertex, and only complete graphs systematically reach this degree. We therefore have a polynomial of degree greater than 1, which corresponds to the general profile observed.

To come back to a complexity expressed only according to the number of vertices, we can consider that $\lvert E\rvert$ is bounded above by $\displaystyle\frac{\lvert V\rvert \times(\lvert V\rvert-1)}{2}$ (which is the case of a complete graph). The asymptotic complexity is therefore $O(\lvert V\rvert^{3})$ (because $\displaystyle\frac{\lvert V\rvert \times(\lvert V\rvert-1)}{2}$ is also bounded by $\lvert V\rvert^{2}$. We find the fact that the graph density influences the computation time of the algorithm.

This approach is very important, the experimental study of the convergence times of an algorithm is absolutely not sufficient, it must be systematically combined with a theoretical study of the temporal complexity.

Do you think we can do better in terms of execution speed ? We could test the adjacency list. What do you think that would produce ?

It has been considered that the temporal complexity depends on the size of the data that is considered. That's true, but a bit imprecise. It depends on the size of the data in the binary coding used to represent the problem. Here, we consider a adjacency matrix (it’s not strictly speaking a binary coding, but it’s obvious that the size in binary will be directly proportional to the size of the matrix), which implies, for example, that the deletion of an edge is done in constant time. It is possible that another representation gives best results (or worse), depending on the operations considered.

To study the algorithm’s behaviour with the representation by adjacency list, it would be necessary to :

1. reimplement the Eulerian cycle algorithm to work with the adjacency list
2. implement a generator creating a graph in the form of adjacency list (or implement a function that translates an adjacency matrix into an adjacency list, in order to be able to use the function that we already have)

We would probably observe a significantly reduced running time, but a similar shape in terms of curve (it is due to the algorithm, not its implementation).

And finally, what do you think we can do to refine our experimental study ?

First of all, as we have noted, the random generation is relatively simplistic, and it could be interesting to improve it (with all the difficulties that have been mentioned above).

Then, it may be interesting to study the impact of certain parameters of the generation of graph on the performance of the algorithm. So far, we have looked at the number of vertices of the graph. But the number of edges is also interesting. On a graph with many edges, the calculation of an Eulerian cycle will inevitably require more computation time. Other parameters could also be studied, such as density (which can also be modified using the probability parameter in the function numpy.random.choice() of the random graph generation), since a very dense graph will also have many edges.

Finally, if the mean time is a good descriptive statistic, it is useful to add the standard deviation to it, as does %%timeit.

4. Conclusion

We have finished this Eulerian cycle ! But we haven’t looked at all the data structures to implement a graph. For example, there is the incidence matrix. Maybe this structure is more efficient to execute our Eulerian cycle algorithm (or less, who knows...).

And if you want to implement the Chinese postman algorithm, you should get it done without a problem. Besides, if you implement an algorithm that applies the Eulerian cycle or the Chinese postman algorithm according to the input data, it could be interesting to study its experimental performance...