<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 2A.algo - Parcourir les rues de Paris
Algorithme de plus courts chemins dans un graphe. Calcul d'un chemin comparé au calcul de tous les chemins les plus courts.
Step1:  Cette idée vient d'une soirée Google Code initiée par Google et à laquelle des élèves de l'ENSAE ont participé. On dispose de la description des rues de Paris (qu'on considèrera comme des lignes droites). On veut déterminer le trajet de huit voitures de telle sorte qu'elles parcourent la ville le plus rapidement possible. On supposera deux cas 
Step2:  Une versions de ce problème est proposée sous forme de challenge 
Step3:  On extrait du fichier l'ensemble des carrefours (vertices) et des rues ou segment de rues (edges).
Step4:  On trace sur un graphique un échantillon des carrefours. On suppose la ville de Paris suffisamment petite et loin des pôles pour considérer les coordonnées comme cartésiennes (et non comme longitude/latitude).
Step5:  Puis on dessine également un échantillon des rues.
Step6:  Petite remarque 
Step7:  Une première solution au premier problème
Ce problème est très similaire à celui du portier chinois. La solution qui suit n'est pas nécessaire la meilleure mais elle donne une idée de ce que peut-être une recherche un peu expérimentale sur le sujet.
Chaque noeud représente un carrefour et chaque rue est un arc reliant des deux carrefours. L'objectif est de parcourir tous les arcs du graphe avec 8 voitures. 
Premiere remarque, l'énoncé ne dit pas qu'il faut parcourir toutes les rues une seule fois. On conçoit aisément que ce serait l'idéal mais on ne sait pas si ce serait possible. Néanmoins, si une telle solution (un chemin passant une et une seule fois par toutes les rues) existe, elle est forcément optimale.
Deuxième remarque, les sens interdits rendent le problème plus complexe. On va dans un premier temps ne pas en tenir compte. On verra comment ajouter la contrainte par la suite et il y a aussi le problème des impasses. On peut néanmoins les contourner en les supprimant du graphe 
Step8:  On ne peut pas trouver une chemin qui parcourt tous les arcs du graphe précédent une et une seule fois. Qu'en est-il du graphe de la ville de Paris ? On compte les noeuds qui ont un nombre pairs et impairs d'arcs les rejoignant (on appelle cela le degré).
Step9:  On remarque que la plupart des carrefours sont le départ de 3 rues. Qu'à cela ne tienne, pourquoi ne pas ajouter des arcs entre des noeuds de degré impair jusqu'à ce qu'il n'y en ait plus que 2. De cette façon, il sera facile de construire un seul chemin parcourant toutes les rues. Comment ajouter ces arcs ? Cela va se faire en deux étapes 
Step10:  Ensuite, on implémente l'algorithme de Bellman-Ford.
Step11:  On s'aperçoit vite que cela va être très très long. On décide alors de ne considérer que les paires de noeuds pour lesquelles la distance à vol d'oiseau est inférieure au plus grand segment de rue ou inférieure à cette distance multipliée par un coefficient.
Step12:  On calcule les arcs admissibles (en espérant que les noeuds de degré impairs seront bien dedans). Cette étape prend quelques minutes 
Step13:  On vérifie que les noeuds de degré impairs font tous partie de l'ensemble des noeuds recevant de nouveaux arcs. La matrice de Bellman envisagera au pire 2.2% de toutes les distances possibles.
Step14:  On continue avec l'algorithme de Bellman-Ford modifié 
Step15:  L'algorithme consiste à regarder les chemins $a \rightarrow b \rightarrow c$ et à comparer s'il est plus rapide que $a \rightarrow c$. 2.6% > 2.2% parce que le filtre est appliqué seulement sur $a \rightarrow b$. Finalement, on considère les arcs ajoutés puis on retire les arcs originaux.
Step16:  Kruskall
On trie les arcs par distance croissante, on enlève les arcs qui ne relient pas des noeuds de degré impair puis on les ajoute un par jusqu'à ce qu'il n'y ait plus d'arc de degré impair.
Step17:  Le nombre de noeuds impairs obtenus à la fin doit être inférieur à 2 pour être sûr de trouver un chemin (mais on choisira 0 pour avoir un circuit eulérien). Mon premier essai n'a pas donné satisfaction (92 noeuds impairs restant) car j'avais choisi un seuil (max_segment / 4) trop petit lors de la sélection des arcs à ajouter. J'ai augmenté le seuil par la suite mais il reste encore 22 noeuds de degré impairs. On a le choix entre augmenter ce seuil mais l'algorithme est déjà long ou chercher dans une autre direction comme laisser l'algorithme de Bellman explorer les noeuds de degré impairs. Ca ne veut pas forcément dire qu'il manque des arcs mais que peut-être ils sont mal choisis. Si l'arc $i \rightarrow j$ est choisi, l'arc $j \rightarrow k$ ne le sera pas car $j$ aura un degré pair. Mais dans ce cas, si l'arc $j \rightarrow k$ était le dernier arc disponible pour combler $k$, on est coincé. On peut augmenter le seuil encore mais cela risquee de prendre du temps et puis cela ne fonctionnerait pas toujours sur tous les jeux de données.
On pourait alors écrire une sorte d'algorithme itératif qui exécute l'algorithme de Bellman, puis lance celui qui ajoute les arcs. Puis on revient au premier en ajoutant plus d'arcs autour des noeuds problèmatique lors de la seconde étape. L'ensemble est un peu long pour tenir dans un notebook mais le code est celui de la fonction eulerien_extension. Je conseille également la lecture de cet article 
Step18:  On enregistre le résultat où on souhaite recommencer à partir de ce moment-là plus tard.
Step19:  Et si vous voulez le retrouver 
Step20:  Chemin Eulérien
A cet instant, on n'a pas vraiment besoin de connaître la longueur du chemin eulérien passant par tous les arcs. Il s'agit de la somme des arcs initiaux et ajoutés (soit environ 334 + 1511). On suppose qu'il n'y qu'une composante connexe. Construire le chemin eulérien fait apparaître quelques difficultés comme la suivante 
Step21:  Quelques algorithmes sont disponibles sur cette page Eulerian_path. L'algorithme de Hierholzer consiste à commencer un chemin eulérien qui peut revenir au premier avant d'avoir tout parcouru. Dans ce cas, on parcourt les noeuds du graphe pour trouver un noeud qui repart ailleurs et qui revient au même noeud. On insert cette boucle dans le chemin initial. Tout d'abord, on construit une structure qui pour chaque noeud associe les noeuds suivant. La fonction euler_path
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
Explanation: 2A.algo - Parcourir les rues de Paris
Algorithme de plus courts chemins dans un graphe. Calcul d'un chemin comparé au calcul de tous les chemins les plus courts.
End of explanation
from jyquickhelper import add_notebook_menu
add_notebook_menu()
Explanation: Cette idée vient d'une soirée Google Code initiée par Google et à laquelle des élèves de l'ENSAE ont participé. On dispose de la description des rues de Paris (qu'on considèrera comme des lignes droites). On veut déterminer le trajet de huit voitures de telle sorte qu'elles parcourent la ville le plus rapidement possible. On supposera deux cas :
Les voitures peuvent être placées n'importe où dans la ville.
Les voitures démarrent et reviennent au même point de départ, le même pour toutes.
Ce notebook décrit comment récupérer les données et propose une solution. Ce problème est plus connu sous le nom du problème du postier chinois ou Route inspection problem pour lequel il existe un algorithme optimal à coût polynomial. Le problème n'est donc pas NP complet.
End of explanation
import pyensae.datasource
data = pyensae.datasource.download_data("paris_54000.zip")
data
Explanation: Une versions de ce problème est proposée sous forme de challenge : City Tour.
Les données
On récupère les données sur Internet.
End of explanation
name = data[0]
with open(name, "r") as f : lines = f.readlines()
    
vertices = []
edges    = [ ]
for i,line in enumerate(lines) :
    spl = line.strip("\n\r").split(" ")
    if len(spl) == 2 :
        vertices.append ( (float(spl[0]), float(spl[1]) ) )
    elif len(spl) == 5 and i > 0:
        v1,v2 = int(spl[0]),int(spl[1])
        ways = int(spl[2]) # dans les deux sens ou pas
        p1 = vertices[v1]
        p2 = vertices[v2]
        edges.append ( (v1,v2,ways,p1,p2) )
    elif i > 0 :
        raise Exception("unable to interpret line {0}: ".format(i) + line)
print("#E=",len(edges), "#V=",len(vertices), ">",max (max( _[0] for _ in edges),  max( _[1] for _ in edges)))
Explanation: On extrait du fichier l'ensemble des carrefours (vertices) et des rues ou segment de rues (edges).
End of explanation
import matplotlib.pyplot as plt
import random
sample = [ vertices[random.randint(0,len(vertices)-1)] for i in range(0,1000)]
plt.plot( [_[0] for _ in sample], [_[1] for _ in sample], ".")
Explanation: On trace sur un graphique un échantillon des carrefours. On suppose la ville de Paris suffisamment petite et loin des pôles pour considérer les coordonnées comme cartésiennes (et non comme longitude/latitude).
End of explanation
sample = [ edges[random.randint(0,len(edges)-1)] for i in range(0,1000)]
for edge in sample:
    plt.plot( [_[0] for _ in edge[-2:]], [_[1] for _ in edge[-2:]], "b-")
Explanation: Puis on dessine également un échantillon des rues.
End of explanation
len ( list(e for e in edges if e[0]==e[1] ))
Explanation: Petite remarque : il n'y a pas de rues reliant le même carrefour :
End of explanation
import networkx as nx
g = nx.Graph()
for i,j in [(1,2),(1,3),(1,4),(2,3),(3,4),(4,5),(5,2),(2,4) ]:
    g.add_edge( i,j )
import matplotlib.pyplot as plt
f, ax = plt.subplots(figsize=(6,3))
nx.draw(g, ax = ax)
Explanation: Une première solution au premier problème
Ce problème est très similaire à celui du portier chinois. La solution qui suit n'est pas nécessaire la meilleure mais elle donne une idée de ce que peut-être une recherche un peu expérimentale sur le sujet.
Chaque noeud représente un carrefour et chaque rue est un arc reliant des deux carrefours. L'objectif est de parcourir tous les arcs du graphe avec 8 voitures. 
Premiere remarque, l'énoncé ne dit pas qu'il faut parcourir toutes les rues une seule fois. On conçoit aisément que ce serait l'idéal mais on ne sait pas si ce serait possible. Néanmoins, si une telle solution (un chemin passant une et une seule fois par toutes les rues) existe, elle est forcément optimale.
Deuxième remarque, les sens interdits rendent le problème plus complexe. On va dans un premier temps ne pas en tenir compte. On verra comment ajouter la contrainte par la suite et il y a aussi le problème des impasses. On peut néanmoins les contourner en les supprimant du graphe : il faut nécessairement faire demi-tour et il n'y a pas de meilleure solution.
Ces deux remarques étant faite, ce problème rappelle un peu le problème des sept ponts de Konigsberg : comment traverser passer par les sept de la ville une et une seule fois. Le mathématicien Euler a répondu à cette question : c'est simple, il suffit que chaque noeud du graphe soit rejoint par un nombre pair de d'arc (= ponts) sauf au plus 2 (les noeuds de départ et d'arrivée). De cette façon, à chaque qu'on rejoint un noeud, il y a toujours une façon de repartir.
End of explanation
nb_edge = { }
for edge in edges :
    v1,v2 = edge[:2]
    nb_edge[v1] = nb_edge.get(v1,0)+1
    nb_edge[v2] = nb_edge.get(v2,0)+1
parite = { }
for k,v in nb_edge.items():
    parite[v] = parite.get(v,0) + 1
[ sorted(parite.items()) ]
Explanation: On ne peut pas trouver une chemin qui parcourt tous les arcs du graphe précédent une et une seule fois. Qu'en est-il du graphe de la ville de Paris ? On compte les noeuds qui ont un nombre pairs et impairs d'arcs les rejoignant (on appelle cela le degré).
End of explanation
def distance(p1,p2):
    return ((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)**0.5
edges = [ edge + (distance( edge[-2],edge[-1]),) for edge in edges]
Explanation: On remarque que la plupart des carrefours sont le départ de 3 rues. Qu'à cela ne tienne, pourquoi ne pas ajouter des arcs entre des noeuds de degré impair jusqu'à ce qu'il n'y en ait plus que 2. De cette façon, il sera facile de construire un seul chemin parcourant toutes les rues. Comment ajouter ces arcs ? Cela va se faire en deux étapes :
On utilise l'algorithme de Bellman-Ford pour construire une matrice des plus courts chemins entre tous les noeuds.
On s'inspire de l'algorithme de poids minimal Kruskal. On trie les arcs par ordre croissant de distance. On ajoute ceux qui réduisent le nombre de noeuds de degré impairs en les prenant dans cet ordre.
Quelques justifications : le meilleur parcours ne pourra pas descendre en deça de la somme des distances des rues puisqu'il faut toutes les parcourir. De plus, s'il existe un chemin qui parcourt toutes les rues, en dédoublant toutes celles parcourues plus d'une fois, il est facile de rendre ce chemin eulérien dans un graphe légèrement modifié par rapport au graphe initial.
Etape 1 : la matrice Bellman
Je ne détaillerai pas trop, la page Wikipedia est assez claire. Dans un premier temps on calcule la longueur de chaque arc (de façon cartésienne). Une autre distance (Haversine) ne changerait pas le raisonnement.
End of explanation
import datetime
init = { (e[0],e[1]) : e[-1] for e in edges }
init.update ( { (e[1],e[0]) : e[-1] for e in edges } )
edges_from = {  }
for e in edges :
    if e[0] not in edges_from : edges_from[e[0]] = []
    if e[1] not in edges_from : edges_from[e[1]] = []
    edges_from[e[0]].append(e)
    edges_from[e[1]].append( (e[1], e[0], e[2], e[4], e[3], e[5] ) )
    
modif = 1
total_possible_edges = len(edges_from)**2
it = 0
while modif > 0 :
    modif = 0
    initc = init.copy()   # to avoid RuntimeError: dictionary changed size during iteration
    s = 0
    for i,d in initc.items() :
        fromi2 = edges_from[i[1]]
        s += d
        for e in fromi2 :
            if i[0] == e[1] : # on fait attention à ne pas ajouter de boucle sur le même noeud
                continue
            new_e = i[0], e[1]
            new_d = d + e[-1]
            if new_e not in init or init[new_e] > new_d :
                init[new_e] = new_d 
                modif += 1
    print(datetime.datetime.now(), "iteration ", it, " modif ", modif, " # ", len(initc),"/",total_possible_edges,"=", 
            "%1.2f" %(len(initc)*100 / total_possible_edges) + "%")
    it += 1
    if it > 6 : break
Explanation: Ensuite, on implémente l'algorithme de Bellman-Ford.
End of explanation
max_segment = max( e[-1] for e in edges )
max_segment
Explanation: On s'aperçoit vite que cela va être très très long. On décide alors de ne considérer que les paires de noeuds pour lesquelles la distance à vol d'oiseau est inférieure au plus grand segment de rue ou inférieure à cette distance multipliée par un coefficient.
End of explanation
possibles = { (e[0],e[1]) : e[-1] for e in edges }
possibles.update ( { (e[1],e[0]) : e[-1] for e in edges } )
initial = possibles.copy()
for i1,v1 in enumerate(vertices) :
    for i2 in range(i1+1,len(vertices)):
        v2 = vertices[i2]
        d = distance(v1,v2)
        if d < max_segment / 2 :         # on ajuste le seuil
            possibles [ i1,i2 ] = d
            possibles [ i2,i1 ] = d
print("original",len(initial),"/",total_possible_edges,"=", len(initial)/total_possible_edges)
print("addition",len(possibles)-len(initial),"/",total_possible_edges,"=", (len(possibles)-len(initial))/total_possible_edges)
Explanation: On calcule les arcs admissibles (en espérant que les noeuds de degré impairs seront bien dedans). Cette étape prend quelques minutes :
End of explanation
allv = { p[0]:True for p in possibles if p not in initial } # possibles est symétrique
for v,p in nb_edge.items() :
    if p % 2 == 1 and v not in allv :
        raise Exception("problème pour le noeud: {0}".format(v))
print("si vous voyez cette ligne, c'est que tout est bon")
Explanation: On vérifie que les noeuds de degré impairs font tous partie de l'ensemble des noeuds recevant de nouveaux arcs. La matrice de Bellman envisagera au pire 2.2% de toutes les distances possibles.
End of explanation
import datetime
init = { (e[0],e[1]) : e[-1] for e in edges }
init.update ( { (e[1],e[0]) : e[-1] for e in edges } )
edges_from = {  }
for e in edges :
    if e[0] not in edges_from : edges_from[e[0]] = []
    if e[1] not in edges_from : edges_from[e[1]] = []
    edges_from[e[0]].append(e)
    edges_from[e[1]].append( (e[1], e[0], e[2], e[4], e[3], e[5] ) )
    
modif = 1
total_possible_edges = len(edges_from)**2
it = 0
while modif > 0 :
    modif = 0
    initc = init.copy()   # to avoid RuntimeError: dictionary changed size during iteration
    s = 0
    for i,d in initc.items() :
        if i not in possibles : 
            continue  # we skip undesired edges ------------------- addition
        fromi2 = edges_from[i[1]]
        s += d
        for e in fromi2 :
            if i[0] == e[1] : # on fait attention à ne pas ajouter de boucle sur le même noeud
                continue
            new_e = i[0], e[1]
            new_d = d + e[-1]
            if new_e not in init or init[new_e] > new_d :
                init[new_e] = new_d 
                modif += 1
    print(datetime.datetime.now(), "iteration ", it, " modif ", modif, " # ", len(initc),"/",total_possible_edges,"=", 
            "%1.2f" %(len(initc)*100 / total_possible_edges) + "%")
    it += 1
    if it > 20 : 
        break
Explanation: On continue avec l'algorithme de Bellman-Ford modifié :
End of explanation
original = { (e[0],e[1]) : e[-1] for e in edges }
original.update ( { (e[1],e[0]) : e[-1] for e in edges } )
additions = { k:v for k,v in init.items() if k not in original }
additions.update( { (k[1],k[0]):v for k,v in additions.items() } )    
Explanation: L'algorithme consiste à regarder les chemins $a \rightarrow b \rightarrow c$ et à comparer s'il est plus rapide que $a \rightarrow c$. 2.6% > 2.2% parce que le filtre est appliqué seulement sur $a \rightarrow b$. Finalement, on considère les arcs ajoutés puis on retire les arcs originaux.
End of explanation
degre = { }
for k,v in original.items() :  # original est symétrique
    degre[k[0]] = degre.get(k[0],0) + 1
tri = [ (v,k) for k,v in additions.items() if degre[k[0]] %2 == 1 and degre[k[1]] %2 == 1  ]
tri.extend( [ (v,k) for k,v in original.items() if degre[k[0]] %2 == 1 and degre[k[1]] %2 == 1  ] )
tri.sort()
impairs = sum ( v%2 for k,v in degre.items()  )
added_edges = [] 
for v,a in tri :
    if degre[a[0]] % 2 == 1 and degre[a[1]] % 2 == 1 :
        # il faut refaire le test car degre peut changer à chaque itération
        degre[a[0]] += 1
        degre[a[1]] += 1
        added_edges.append ( a + (v,) )
        impairs -= 2
        if impairs <= 0 :
            break
# on vérifie
print("nb degré impairs",impairs, "nombre d'arcs ajoutés",len(added_edges))
print("longueur ajoutée ", sum( v for a,b,v in added_edges ))
print("longueur initiale ", sum( e[-1] for e in edges ))
Explanation: Kruskall
On trie les arcs par distance croissante, on enlève les arcs qui ne relient pas des noeuds de degré impair puis on les ajoute un par jusqu'à ce qu'il n'y ait plus d'arc de degré impair.
End of explanation
from ensae_teaching_cs.special.rues_paris import eulerien_extension, distance_paris,get_data
print("data")
edges = get_data()
print("start, nb edges", len(edges))
added = eulerien_extension(edges, distance=distance_paris)
print("end, nb added", len(added))
Explanation: Le nombre de noeuds impairs obtenus à la fin doit être inférieur à 2 pour être sûr de trouver un chemin (mais on choisira 0 pour avoir un circuit eulérien). Mon premier essai n'a pas donné satisfaction (92 noeuds impairs restant) car j'avais choisi un seuil (max_segment / 4) trop petit lors de la sélection des arcs à ajouter. J'ai augmenté le seuil par la suite mais il reste encore 22 noeuds de degré impairs. On a le choix entre augmenter ce seuil mais l'algorithme est déjà long ou chercher dans une autre direction comme laisser l'algorithme de Bellman explorer les noeuds de degré impairs. Ca ne veut pas forcément dire qu'il manque des arcs mais que peut-être ils sont mal choisis. Si l'arc $i \rightarrow j$ est choisi, l'arc $j \rightarrow k$ ne le sera pas car $j$ aura un degré pair. Mais dans ce cas, si l'arc $j \rightarrow k$ était le dernier arc disponible pour combler $k$, on est coincé. On peut augmenter le seuil encore mais cela risquee de prendre du temps et puis cela ne fonctionnerait pas toujours sur tous les jeux de données.
On pourait alors écrire une sorte d'algorithme itératif qui exécute l'algorithme de Bellman, puis lance celui qui ajoute les arcs. Puis on revient au premier en ajoutant plus d'arcs autour des noeuds problèmatique lors de la seconde étape. L'ensemble est un peu long pour tenir dans un notebook mais le code est celui de la fonction eulerien_extension. Je conseille également la lecture de cet article : Efficient Algorithms for Eulerian Extension (voir également On Making Directed Graphs Eulerian). L'exécution qui suit prend une vingtaine de minutes.
End of explanation
with open("added.txt","w") as f : f.write(str(added))
Explanation: On enregistre le résultat où on souhaite recommencer à partir de ce moment-là plus tard.
End of explanation
from ensae_teaching_cs.data import added
data = added()
from ensae_teaching_cs.special.rues_paris import eulerien_extension, distance_paris, get_data
edges = get_data()
with open(data, "r") as f: text = f.read()
added_edges = eval(text)
Explanation: Et si vous voulez le retrouver :
End of explanation
from pyquickhelper.helpgen import NbImage
NbImage("euler.png")
Explanation: Chemin Eulérien
A cet instant, on n'a pas vraiment besoin de connaître la longueur du chemin eulérien passant par tous les arcs. Il s'agit de la somme des arcs initiaux et ajoutés (soit environ 334 + 1511). On suppose qu'il n'y qu'une composante connexe. Construire le chemin eulérien fait apparaître quelques difficultés comme la suivante : on parcourt le graphe dans un sens mais on peut laisser de côté une partie du chemin et créer une seconde composante connexe.
End of explanation
from ensae_teaching_cs.special.rues_paris import euler_path
path = euler_path(edges, added_edges)
path[:5]
Explanation: Quelques algorithmes sont disponibles sur cette page Eulerian_path. L'algorithme de Hierholzer consiste à commencer un chemin eulérien qui peut revenir au premier avant d'avoir tout parcouru. Dans ce cas, on parcourt les noeuds du graphe pour trouver un noeud qui repart ailleurs et qui revient au même noeud. On insert cette boucle dans le chemin initial. Tout d'abord, on construit une structure qui pour chaque noeud associe les noeuds suivant. La fonction euler_path
End of explanation
<END_TASK>