Coverage for src/ensae_teaching_cs/special/tsp_kruskal.py : 93%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2"""
3@file
4@brief Implémente un algorithme qui cherche le plus court chemin passant
5par tous les noeuds d'un graphe (TSP). Applique un algorithme de Kruskal
6puis cherche à améliorer le chemin localement.
7Voir :ref:`l-tsp_kruskal`. La fonction principale est
8@see fn tsp_kruskal_algorithm.
9"""
10import functools
11import random
12import math
13import os
14from pyquickhelper.loghelper import noLOG
15from .tsp_bresenham import draw_line
16from ..helpers.pygame_helper import wait_event, empty_main_loop
19def construit_ville(n, x=1000, y=700):
20 """
21 Tire aléatoirement *n* villes dans un carrée ``x * y``, on choisit
22 ces villes de sortent qu'elles ne soient pas trop proches."""
23 # deux villes ne pourront pas être plus proches que mind
24 mind = math.sqrt(x * x + y * y) / (n * 0.75)
25 # liste vide
26 lv = []
27 while n > 0:
28 # on tire aléatoirement les coordonnées d'une ville
29 xx = x * random.random()
30 yy = y * random.random()
31 # on vérifie qu'elle n'est pas trop proche d'aucune autre ville
32 ajout = True
33 for t in lv:
34 d1 = t[0] - xx
35 d2 = t[1] - yy
36 d = math.sqrt(d1 * d1 + d2 * d2)
37 if d < mind:
38 ajout = False # ville trop proche
39 # si la ville n'est pas trop proche des autres, on l'ajoute à la liste
40 if ajout:
41 lv.append((xx, yy))
42 n -= 1 # une ville en moins à choisir
43 return lv
46def distance_euclidian(p1, p2):
47 """
48 Calcule la distance entre deux villes.
49 """
50 x = p1[0] - p2[0]
51 y = p1[1] - p2[1]
52 return math.sqrt(x * x + y * y)
55def vecteur_points(p1, p2):
56 """
57 Retourne le vecteur entre les points *p1* et *p2*.
58 """
59 return (p2[0] - p1[0], p2[1] - p1[1])
62def vecteur_norme(vec):
63 """
64 Retourne la norme d'un vecteur.
65 """
66 return math.sqrt(vec[0] * vec[0] + vec[1] * vec[1])
69def vecteur_cosinus(vec1, vec2):
70 """
71 Retourne le cosinus entre deux vecteurs,
72 utilise le produit scalaire.
73 """
74 norm1 = vecteur_norme(vec1)
75 norm2 = vecteur_norme(vec2)
76 if norm1 == 0:
77 return 1
78 if norm2 == 0:
79 return 1
80 scal = vec1[0] * vec2[0] + vec1[1] * vec2[1]
81 return scal / (norm1 * norm2)
84def vecteur_sinus(vec1, vec2):
85 """
86 Retourne le sinus entre deux vecteurs,
87 utilise le produit vectoriel.
88 """
89 norm1 = vecteur_norme(vec1)
90 norm2 = vecteur_norme(vec2)
91 if norm1 == 0:
92 return 0
93 if norm2 == 0:
94 return 0
95 scal = vec1[0] * vec2[1] - vec1[1] * vec2[0]
96 return scal / (norm1 * norm2)
99def oppose_vecteur(vec):
100 """
101 retourne le vecteur opposé.
102 """
103 return (-vec[0], -vec[1])
106def repartition_zone(villes, zone_taille, ask_zone=False):
107 """
108 Répartit les villes en zones, retourne les villes rangées par zones,
109 chaque éléments zones [z][k] contient :
111 - les coordonnées de la ville
112 - ses coordonnées en zone, (zx, zy)
113 - son indice dans la liste villes
115 La fonction retourne également le nombre de zones
116 selon l'axe des abscisses et l'axe des ordonnées,
117 retourne aussi le nombre de zones, si *ask_zone* est True,
118 retourne un paramètre supplémentaire : *zone*.
119 """
120 mx = min(v[0] for v in villes)
121 my = min(v[1] for v in villes)
122 X = max((v[0] - mx) / zone_taille for v in villes)
123 Y = max((v[1] - my) / zone_taille for v in villes)
124 X = int(X) + 1
125 Y = int(Y) + 1
127 # attribution des zones
128 zone = []
129 Zmax = 0
130 for i in range(0, len(villes)):
131 v = villes[i]
132 x = int((v[0] - mx) / zone_taille)
133 y = int((v[1] - my) / zone_taille)
134 z = int(y * X + x)
135 Zmax = max(z, Zmax)
136 zone.append((z, v, (x, y), i))
138 # rangement par zone
139 Zmax += 1
140 zones = [[] for i in range(0, Zmax)]
141 for z in zone:
142 zones[z[0]].append((z[1], z[2], z[3]))
144 if ask_zone:
145 return zones, X, Y, mx, my, Zmax, zone
146 return zones, X, Y, mx, my, Zmax
149def voisinage_zone(z, Zmax, X, Y):
150 """
151 Retourne la liste des voisins d'une zone *z*
152 sachant qu'il y a *X* zones sur l'axe des abscisses
153 et *Y* zones sur l'axe des ordonnées,
154 *Zmax* est le nombre de zones,
155 inclus *z* dans cette liste.
156 """
157 x = z % X
158 y = z // X
159 voisin_ = [z]
160 if x > 0:
161 voisin_.append(z - 1)
162 if x < X:
163 voisin_.append(z + 1)
164 if y > 0:
165 voisin_.append(z - X)
166 if y < Y:
167 voisin_.append(z + X)
168 if x > 0 and y > 0:
169 voisin_.append(z - 1 - X)
170 if x > 0 and y < Y:
171 voisin_.append(z - 1 + X)
172 if x < X and y > 0:
173 voisin_.append(z + 1 - X)
174 if x < X and y < Y:
175 voisin_.append(z + 1 + X)
176 voisin = [int(i) for i in voisin_ if Zmax > i >= 0]
177 return voisin
180def arbre_poids_minimal(villes, zone_taille, distance):
181 """
182 Construit l'arbre de poids minimal, retourne une liste de
183 listes, chaque sous-liste associée à une ville contient la liste des ses voisins,
184 *zone_taille* permet de découper l'image en zones,
185 les distances ne seront calculées que si
186 deux éléments sont dans la même zone ou dans une zone voisine.
188 @param villes list of tuples (tuple = coordinates)
189 @param zone_taille @see fn repartition_zone
190 @param distance distance function which returns the distance between two
191 elements
192 @return list of lists: each sublist *r[i]* contains the indexes of
193 neighbors of node *i* so that the whole graph is
194 only one connected component
195 """
197 def tri_distance(u, v):
198 if u[2] < v[2]:
199 return -1
200 elif u[2] > v[2]:
201 return 1
202 else:
203 return 0
205 rz = repartition_zone(villes, zone_taille=zone_taille)
206 zones, X, Y, mx, my, Zmax = rz[:6]
208 # calcul des distances
209 li = []
210 for z in range(0, len(zones)):
211 voisin = voisinage_zone(z, Zmax, X, Y)
212 for v in zones[z]:
213 for zz in voisin:
214 for u in zones[zz]:
215 d = distance(v[0], u[0])
216 li.append((v[2], u[2], d))
218 # tri
219 li = list(sorted(li, key=functools.cmp_to_key(tri_distance)))
221 # nombre de composantes connexes
222 nb_comp = len(villes)
224 # indice de la composante d'une ville
225 num_comp = [i for i in range(0, len(villes))]
227 # liste des voisins pour chaque ville
228 arbre = [[] for i in range(0, len(villes))]
230 # liste des villes par composante connexe
231 list_comp = [[i] for i in range(0, len(villes))]
233 while nb_comp > 1:
234 iii = 0
235 for c in li:
236 iii += 1
237 i, j = c[0], c[1]
238 if num_comp[i] != num_comp[j]:
239 # on relie les villes i et j car elles appartiennent
240 # à des composantes connexes différentes
241 arbre[i].append(j) # i est voisine de j
242 arbre[j].append(i) # j est voisine de i
243 cl = num_comp[i] # composante connexe restante
244 # composante connexe à agréger à la précédente
245 ki = num_comp[j]
246 for k in list_comp[ki]:
247 num_comp[k] = cl
248 list_comp[cl].append(k)
249 list_comp[ki] = []
250 nb_comp -= 1 # une composante connexe en moins
252 if nb_comp <= 1:
253 break # il n'y a plus qu'une seule composante connexe, inutile de continuer
255 if nb_comp > 1:
256 # it usually means that zone_taille is too small and some edges
257 # we find lost connected components
258 # so for these, assuming they are not too many
259 # we look for the closest point outside the connected component
260 first_count = min((len(l), i)
261 for i, l in enumerate(list_comp) if len(l) > 0)
262 comp = first_count[1]
263 city = list_comp[comp][random.randint(0, len(list_comp[comp]) - 1)]
264 # city is not the best choice, just a random one
265 dist = min((distance(villes[city], v), i) for i, v in enumerate(villes)
266 if city != i and num_comp[i] != num_comp[city])
267 li = [(city, dist[1])]
269 return dict(arbre=arbre, X=X, Y=Y, mx=mx, my=my, Zmax=Zmax)
272def circuit_eulerien(villes, arbre, screen, pygame, fLOG):
273 """
274 Définit un circuit eulérien, villes contient la liste des villes,
275 tandis que arbre est une liste de listes, ``arbre[i]`` est la liste des villes
276 connectées à la ville *i*,
277 on suppose que arbre est un graphe de poids minimal non orienté,
278 l'algorithme ne marche pas s'il existe des villes confondues,
279 un circuit eulérien passe par tous les arêtes une et une seule fois.
280 """
282 # on choisit une ville qui est une extrémité et parmi celle-là on la
283 # choisit au hasard
284 has = []
285 for i in range(0, len(villes)):
286 n = len(arbre[i])
287 if n == 1:
288 has.append(i)
290 bm = random.randint(0, len(has) - 1)
291 bm = has[bm]
293 # vecteur, le circuit eulérien contient
294 # nécessairement 2 * len (villes) noeuds puisque c'est
295 # le graphe eulérien d'un arbre de poids minimal non orienté
296 vec = (1, 1)
297 chemin = [bm]
298 done = set()
299 done.add(bm)
300 iter = []
301 while len(done) < len(villes):
302 iter.append(len(done))
303 if len(iter) % 1000 == 0:
304 fLOG(" circuit_eulerien: iter={0} len(done)={1} len(villes)={2}".format(
305 len(iter), len(done), len(villes)))
306 if len(done) == iter[-1000]:
307 # there is apparently something wrong
308 break
309 v = villes[bm]
310 ma = - math.pi - 1
311 bvec = vec
312 opvec = oppose_vecteur(vec)
313 bl = None
314 for k in range(0, len(arbre[bm])):
315 la = arbre[bm][k]
316 vec2 = vecteur_points(v, villes[la])
317 if vec2 == (0.0, 0.0):
318 # same point, we keep the same direction
319 if la not in done:
320 bl = la
321 bvec = vec2
322 # no need to go further if the points are equal
323 break
324 # we skip
325 continue
326 if opvec == vec2:
327 angle = -math.pi
328 else:
329 cos = vecteur_cosinus(vec, vec2)
330 sin = vecteur_sinus(vec, vec2)
331 angle = math.atan2(sin, cos)
332 if angle > ma:
333 ma = angle
334 bl = la
335 bvec = vec2
337 if bl is not None:
338 if bl not in done:
339 chemin.append(bl)
340 done.add(bl)
341 bm = bl
342 if bvec != (0.0, 0.0):
343 vec = bvec
344 else:
345 # something is wrong (it might an issue with duplicated points)
346 rows = []
347 for i, p in enumerate(villes):
348 rows.append("p{0}: {1},{2}".format(i, p[0], p[1]))
349 for i, c in enumerate(chemin):
350 rows.append("c{0}: i={1} -> {2},{3}".format(i,
351 c, villes[c][0], villes[c][1]))
352 rows.append("bm={0} ma={1} bvec={2} vec={3} bl={4}".format(
353 bm, ma, vec2, vec, bl))
354 rows.append("arbre[{0}]={1}".format(bm, arbre[bm]))
355 rows.append("arbre[{0}]={1}".format(
356 arbre[bm][0], arbre[arbre[bm][0]]))
357 mes = "\n".join(rows)
358 raise Exception("this case should not happen\n" + mes)
360 if len(done) < len(villes):
361 # something is wrong (it might an issue with duplicated points)
362 rows = []
363 for i, p in enumerate(villes):
364 rows.append("p{0}: {1},{2}".format(i, p[0], p[1]))
365 for i, c in enumerate(chemin):
366 rows.append("c{0}: i={1} -> {2},{3}".format(i,
367 c, villes[c][0], villes[c][1]))
368 rows.append("bm={0} ma={1} bvec={2} vec={3} bl={4}".format(
369 bm, ma, vec2, vec, bl))
370 mes = "\n".join(rows)
371 raise Exception("circuit_eulerien cannot give a path:\n" + mes)
373 return chemin
376def circuit_hamiltonien(chemin):
377 """
378 Extrait un circuit hamiltonien depuis un circuit eulérien,
379 passe par tous les sommets une et une seule fois.
380 """
381 nb = max(chemin) + 1
382 res = []
383 coche = [False for i in range(0, nb)]
384 for c in chemin:
385 if coche[c]:
386 continue
387 res.append(c)
388 coche[c] = True
389 return res
392def equation_droite(p1, p2):
393 """
394 retourne l'équation d'une droite passant par p1 et p2,
395 ax + by + c = 0, retourne les coefficients a,b,c
396 """
397 vec = vecteur_points(p1, p2)
398 a = vec[1]
399 b = - vec[0]
400 c = - a * p1[0] - b * p1[1]
401 return a, b, c
404def evaluation_droite(a, b, c, p):
405 """
406 L'équation d'une droite est : ``ax + by + c``, retourne la valeur
407 de cette expression au point *p*.
408 """
409 return a * p[0] + b * p[1] + c
412def intersection_segment(p1, p2, p3, p4):
413 """
414 Dit si les segments *[p1 p2]* et *[p3 p4]* ont une intersection,
415 ne retourne pas l'intersection.
416 """
417 # équation de la droite (p1 p2)
418 a1, b1, c1 = equation_droite(p1, p2)
419 a2, b2, c2 = equation_droite(p3, p4)
420 s1 = evaluation_droite(a2, b2, c2, p1)
421 s2 = evaluation_droite(a2, b2, c2, p2)
422 s3 = evaluation_droite(a1, b1, c1, p3)
423 s4 = evaluation_droite(a1, b1, c1, p4)
424 return s1 * s2 <= 0 and s3 * s4 <= 0
427def longueur_chemin(chemin, distance):
428 """
429 Retourne la longueur d'un chemin.
430 """
431 s = 0
432 nb = len(chemin)
433 for i in range(0, nb):
434 s += distance(chemin[i], chemin[(i + 1) % nb])
435 return s
438def retournement_essai(chemin, i, j, negligeable=1e-5, distance=None):
439 """
440 Dit s'il est judicieux de parcourir le chemin entre les sommets *i* et *j*
441 en sens inverse, si c'est judicieux, change le chemin et retourne True,
442 sinon, retourne False,
443 si *j < i*, alors la partie à inverser est
444 *i*, *i+1*, *i+2*, ..., *len(chemin)*,
445 *-1*, *0*, *1*, ..., *j*.
446 """
448 nb = len(chemin)
449 if i == j:
450 return False
451 if j - i == -1:
452 return False
453 if j - i - nb == -1:
454 return False
456 ia = (i - 1 + nb) % nb
457 ja = (j + 1) % nb
458 # arcs enlevés
459 d_ia_i = distance(chemin[ia], chemin[i])
460 d_j_ja = distance(chemin[j], chemin[ja])
461 # arcs ajoutés
462 d_ia_j = distance(chemin[ia], chemin[j])
463 d_i_ja = distance(chemin[i], chemin[ja])
464 # amélioration ?
465 d = d_ia_j + d_i_ja - d_j_ja - d_ia_i
466 if d >= -negligeable:
467 return False
469 # si amélioration, il faut retourner le chemin entre les indices i et j
470 jp = j
471 if jp < i:
472 jp = j + nb
473 ip = i
475 while ip < jp:
476 i = ip % nb
477 j = jp % nb
478 ech = chemin[i]
479 chemin[i] = chemin[j]
480 chemin[j] = ech
481 ip = ip + 1
482 jp = jp - 1
484 return True
487def retournement(chemin, taille, fLOG, distance):
488 """
489 Amélioration du chemin par un algorithme simple,
490 utilise des retournements de taille au plus <taille>,
491 retourne le nombre de modifications.
492 """
494 # traitement des petits retournements
495 nb = len(chemin)
496 nb_change = 1
497 nbtout = 0
498 retour = {}
499 while nb_change > 0:
500 nb_change = 0
501 for t in range(1, taille + 1):
502 retour[t] = 0
503 for i in range(0, nb):
504 j = (i + t) % nb
505 b = retournement_essai(chemin, i, j, distance=distance)
506 if b:
507 retour[t] += 1
508 nb_change += 1
509 nbtout += nb_change
510 fLOG("nombre de retournements %d longueur : \t %10.0f --- \t"
511 % (nbtout, longueur_chemin(chemin, distance)), " --- : ", retour)
512 return nbtout
515def echange_position_essai(chemin, a, b, x, inversion, negligeable=1e-5, distance=None):
516 """
517 Echange la place des villes ka et kb pour les placer entre les villes *i* et *i+1*,
518 si inversion est True, on inverse également le chemin inséré, si inversion est False,
519 on ne l'inverse pas,
520 si cela améliore, déplace les villes et retourne True, sinon, retourne False.
521 """
523 nb = len(chemin)
524 xa = (x + 1) % nb
525 ka = (a - 1 + nb) % nb
526 kb = (b + 1) % nb
528 if not inversion:
530 if x == ka:
531 return False
532 if x == kb:
533 return False
534 if xa == ka:
535 return False
536 if b < a:
537 if a <= x <= b + nb:
538 return False
539 elif a <= x <= b:
540 return False
541 if b < a:
542 if a <= x + nb <= b + nb:
543 return False
544 elif a <= x + nb <= b:
545 return False
547 # arcs enlevés
548 d_x_xa = distance(chemin[x], chemin[xa])
549 d_ka_a = distance(chemin[ka], chemin[a])
550 d_b_kb = distance(chemin[b], chemin[kb])
551 # arcs ajoutés
552 d_ka_kb = distance(chemin[ka], chemin[kb])
553 d_x_a = distance(chemin[x], chemin[a])
554 d_b_xa = distance(chemin[b], chemin[xa])
556 d = d_ka_kb + d_x_a + d_b_xa - d_x_xa - d_ka_a - d_b_kb
557 if d >= -negligeable:
558 return False
560 # villes à déplacer
561 ech = []
562 bp = b
563 if bp < a:
564 bp = b + nb
565 for i in range(a, bp + 1):
566 ech.append(chemin[i % nb])
567 diff = bp - a + 1
569 xp = x
570 if xp < b:
571 xp += nb
573 for le in range(b + 1, xp + 1):
574 ll = le % nb
575 bp = (a + le - b - 1) % nb
576 chemin[bp] = chemin[ll]
578 for le in range(0, len(ech)):
579 chemin[(x + le - diff + 1 + nb) % nb] = ech[le]
581 return True
583 else:
585 if x == ka:
586 return False
587 if x == kb:
588 return False
589 if xa == ka:
590 return False
591 if b < a:
592 if a <= x <= b + nb:
593 return False
594 elif a <= x <= b:
595 return False
596 if b < a:
597 if a <= x + nb <= b + nb:
598 return False
599 elif a <= x + nb <= b:
600 return False
602 # arcs enlevés
603 d_x_xa = distance(chemin[x], chemin[xa])
604 d_ka_a = distance(chemin[ka], chemin[a])
605 d_b_kb = distance(chemin[b], chemin[kb])
606 # arcs ajoutés
607 d_ka_kb = distance(chemin[ka], chemin[kb])
608 d_x_b = distance(chemin[x], chemin[b])
609 d_a_xa = distance(chemin[a], chemin[xa])
611 d = d_ka_kb + d_x_b + d_a_xa - d_x_xa - d_ka_a - d_b_kb
612 if d >= -negligeable:
613 return False
615 # villes à déplacer
616 ech = []
617 bp = b
618 if bp < a:
619 bp = b + nb
620 for i in range(a, bp + 1):
621 ech.append(chemin[i % nb])
622 ech.reverse()
623 diff = bp - a + 1
625 xp = x
626 if xp < b:
627 xp += nb
629 for le in range(b + 1, xp + 1):
630 ll = le % nb
631 bp = (a + le - b - 1) % nb
632 chemin[bp] = chemin[ll]
634 for le in range(0, len(ech)):
635 chemin[(x + le - diff + 1 + nb) % nb] = ech[le]
637 return True
640def dessin_arete_zone(chemin, taille_zone, X, Y):
641 """
642 Retourne une liste de listes de listes,
643 ``res[i][j]`` est une liste des arêtes passant près de la zone ``(x,y) = [i][j]``,
644 si *k* in ``res[i][j]``, alors l'arête *k*, *k+1* est dans la zone *(i,j)*,
645 *X* est le nombre de zones horizontalement, *Y* est le nombre de zones verticalement,
646 *taille_zone* est la longueur du côté du carré d'une zone.
647 """
648 res = [[[] for j in range(0, Y + 1)] for i in range(0, X + 1)]
649 nb = len(chemin)
650 mx = min(_[0] for _ in chemin)
651 my = min(_[1] for _ in chemin)
652 for i in range(0, nb):
653 a = chemin[i]
654 b = chemin[(i + 1) % nb]
655 x1, x2 = int(
656 (a[0] - mx) // taille_zone), int((b[0] - mx) // taille_zone)
657 y1, y2 = int(
658 (a[1] - my) // taille_zone), int((b[1] - my) // taille_zone)
659 line = draw_line(x1, y1, x2, y2)
660 for x, y in line:
661 res[x][y].append(i)
662 return res
665def voisinage_zone_xy(x, y, X, Y):
666 """
667 Retourne la liste des voisins d'une zone *(x,y)*
668 sachant qu'il y a *X* zones sur l'axe des abscisses
669 et *Y* zones sur l'axe des ordonnées,
670 inclus *z* dans cette liste
671 """
672 voisin = [(x, y)]
673 if x > 0:
674 voisin.append((x - 1, y))
675 if x < X - 1:
676 voisin.append((x + 1, y))
677 if y > 0:
678 voisin.append((x, y - 1))
679 if y < Y - 1:
680 voisin.append((x, y + 1))
681 if x > 0 and y > 0:
682 voisin.append((x - 1, y - 1))
683 if x > 0 and y < Y - 1:
684 voisin.append((x - 1, y + 1))
685 if x < X - 1 and y > 0:
686 voisin.append((x + 1, y - 1))
687 if x < X - 1 and y < Y - 1:
688 voisin.append((x + 1, y + 1))
689 return voisin
692def echange_position(chemin, taille, taille_zone, X, Y, grande=0.5, fLOG=None, distance=None):
693 """
694 Regarde si on ne peut pas déplacer un segment de longueur taille
695 pour supprimer les arêtes les plus longues,
696 au maximum <grande> longues arêtes,
697 retourne le nombre de changement effectués,
698 *X* est le nombre de zones horizontalement,
699 *Y* est le nombre de zones verticalement,
700 *taille_zone* est la longueur d'un côté du carré d'une zone.
701 """
703 nb = len(chemin)
705 def tri_arete(x, y):
706 """pour trier la liste l par ordre décroissant"""
707 if x[2] < y[2]:
708 return 1
709 elif x[2] > y[2]:
710 return -1
711 else:
712 return 0
714 tmx = min(v[0] for v in chemin)
715 tmy = min(v[1] for v in chemin)
717 # list des arêtes triés par ordre décroissant
718 la = []
719 for i in range(0, nb):
720 im = (i + 1) % nb
721 la.append((i, im, distance(chemin[i], chemin[im])))
722 la = list(sorted(la, key=functools.cmp_to_key(tri_arete)))
724 # zone associée à chaque arête
725 zone = dessin_arete_zone(chemin, taille_zone, X, Y)
727 dseuil = la[int(nb * grande)][2]
728 nbtout = 0
729 nb_change = 0
730 iarete = 0
731 retour = {}
732 for t in range(1, taille + 1):
733 retour[t] = 0
735 while iarete < nb:
736 nb_change = 0
737 arete = la[iarete]
738 iarete += 1
739 x = arete[0]
740 xm = arete[1]
741 a = chemin[x]
742 b = chemin[xm]
743 d = distance(a, b)
744 if d < dseuil:
745 break # arête trop petite
747 # zone traversée par la ligne
748 x1, x2 = (int((a[0] - tmx) // taille_zone),
749 int((b[0] - tmx) // taille_zone))
750 y1, y2 = (int((a[1] - tmy) // taille_zone),
751 int((b[1] - tmy) // taille_zone))
752 ens = draw_line(x1, y1, x2, y2)
753 ville = []
754 for k, l in ens:
755 voisin = voisinage_zone_xy(k, l, X, Y)
756 for u, v in voisin:
757 ville.extend(zone[u][v])
759 # on supprime les doubles
760 ville.sort()
761 if len(ville) == 0:
762 continue
763 sup = []
764 mx = -1
765 for v in ville:
766 if v == mx:
767 sup.append(v)
768 mx = v
769 for s in sup:
770 ville.remove(s)
772 # on étudie les possibilités de casser l'arête (x,xm) aux alentours des villes
773 # comprise dans l'ensemble ville
774 for t in range(1, taille + 1):
776 for i in ville:
778 # on essaye d'insérer le sous-chemin (x- t + 1 + nb) --> x
779 # au milieu de l'arête i,i+1
780 b = echange_position_essai(
781 chemin, (x - t + 1 + nb) % nb, x, i, False, distance=distance)
782 if b:
783 nb_change += 1
784 retour[t] += 1
785 continue
787 # on essaye d'insérer le sous-chemin (xm+ t - 1) --> xm
788 # au milieu de l'arête i,i+1
789 b = echange_position_essai(
790 chemin, (xm + t - 1) % nb, xm, i, False, distance=distance)
791 if b:
792 nb_change += 1
793 retour[t] += 1
794 continue
796 # on essaye de casser l'arête x,xm en insérant
797 # le sous-chemin i --> (i+t) % nb
798 b = echange_position_essai(
799 chemin, i, (i + t) % nb, x, False, distance=distance)
800 if b:
801 nb_change += 1
802 retour[t] += 1
803 continue
804 # idem
805 b = echange_position_essai(
806 chemin, i, (i + t) % nb, x, True, distance=distance)
807 if b:
808 retour[t] += 1
809 nb_change += 1
810 continue
811 # idem
812 b = echange_position_essai(
813 chemin, (i - t + nb) % nb, i, x, False, distance=distance)
814 if b:
815 nb_change += 1
816 retour[t] += 1
817 continue
818 # idem
819 b = echange_position_essai(
820 chemin, (i - t + nb) % nb, i, x, True, distance=distance)
821 if b:
822 retour[t] += 1
823 nb_change += 1
824 continue
826 nbtout += nb_change
828 fLOG("nombre de déplacements %d longueur : \t %10.0f --- \t"
829 % (nbtout, longueur_chemin(chemin, distance=distance)), " --- : ", retour)
830 return nbtout
833def supprime_croisement(chemin, taille_zone, X, Y, fLOG, distance=None):
834 """
835 Supprime les croisements d'arêtes,
836 retourne le nombre de changement effectués,
837 *X* est le nombre de zones horizontalement,
838 *Y* est le nombre de zones verticalement,
839 *taille_zone* est la longueur d'un côté du carré d'une zone
840 """
842 nb = len(chemin)
843 tmx = min(v[0] for v in chemin)
844 tmy = min(v[1] for v in chemin)
846 # zone associée à chaque arête
847 zone = dessin_arete_zone(chemin, taille_zone, X, Y)
848 nbtout = 0
850 for i in range(0, nb):
851 im = (i + 1) % nb
852 a = chemin[i]
853 b = chemin[im]
855 # zone traversée par la ligne
856 x1, x2 = (int((a[0] - tmx) // taille_zone),
857 int((b[0] - tmx) // taille_zone))
858 y1, y2 = (int((a[1] - tmy) // taille_zone),
859 int((b[1] - tmy) // taille_zone))
860 ens = draw_line(x1, y1, x2, y2)
861 ville = []
862 for k, l in ens:
863 voisin = voisinage_zone_xy(k, l, X, Y)
864 for u, v in voisin:
865 ville.extend(zone[u][v])
867 # on supprime les doubles
868 ville.sort()
869 if len(ville) == 0:
870 continue
871 sup = []
872 mx = -1
873 for v in ville:
874 if v == mx:
875 sup.append(v)
876 mx = v
877 for s in sup:
878 ville.remove(s)
880 nb_change = 0
881 for v in ville:
882 b = retournement_essai(chemin, i, v, distance=distance)
883 if b:
884 nb_change += 1
885 continue
886 b = retournement_essai(chemin, im, v, distance=distance)
887 if b:
888 nb_change += 1
889 continue
891 nbtout += nb_change
893 fLOG("nombre de croisements %d longueur : \t %10.0f --- \t"
894 % (nbtout, longueur_chemin(chemin, distance=distance)))
895 return nbtout
898def amelioration_chemin(chemin, taille_zone, X, Y, taille=10, screen=None,
899 fLOG=None, pygame=None, max_iter=None, images=None,
900 distance=None):
901 """
902 Amélioration du chemin par un algorithme simple,
903 utilise des retournements de taille au plus *taille*,
904 traite les arcs qui se croisent,
905 traite les grands arcs, utilise un quadrillage de taille *window*,
906 *X* est le nombre de zones horizontalement,
907 *Y* est le nombre de zones verticalement,
908 *taille_zone* est la longueur d'un côté du carré d'une zone.
909 """
911 white = 255, 255, 255
913 if pygame is not None and images is not None:
914 images.append(screen.copy())
916 # première étape rapide
917 iter = 0
918 nb = 1
919 while nb > 0 and (max_iter is None or iter < max_iter):
920 nb = retournement(chemin, taille, fLOG=fLOG, distance=distance)
921 if screen is not None:
922 screen.fill(white)
923 display_chemin(chemin, 0, screen, pygame=pygame)
924 pygame.display.flip()
925 if images is not None:
926 images.append(screen.copy())
927 empty_main_loop(pygame)
928 iter += 1
930 # amélioration
931 nb = 1
932 while nb > 0 and (max_iter is None or iter < max_iter):
933 nb = retournement(chemin, taille, fLOG=fLOG, distance=distance)
934 if screen is not None:
935 screen.fill(white)
936 display_chemin(chemin, 0, screen=screen, pygame=pygame)
937 pygame.display.flip()
938 if images is not None:
939 images.append(screen.copy())
940 empty_main_loop(pygame)
941 nb += echange_position(chemin, taille // 2,
942 taille_zone, X, Y, fLOG=fLOG,
943 distance=distance)
944 if screen is not None:
945 screen.fill(white)
946 display_chemin(chemin, 0, screen=screen, pygame=pygame)
947 pygame.display.flip()
948 if images is not None:
949 images.append(screen.copy())
950 empty_main_loop(pygame)
951 nb += supprime_croisement(chemin, taille_zone,
952 X, Y, fLOG=fLOG, distance=distance)
953 if screen is not None:
954 screen.fill(white)
955 display_chemin(chemin, 0, screen=screen, pygame=pygame)
956 pygame.display.flip()
957 if images is not None:
958 images.append(screen.copy())
959 empty_main_loop(pygame)
960 iter += 1
963def tsp_kruskal_algorithm(points, size=20, length=10, max_iter=None,
964 fLOG=noLOG, pygame=None, screen=None, images=None,
965 distance=None):
966 """
967 Finds the shortest path going through all points,
968 points require to be a 2 dimensional space.
970 @param points list 2-tuple (X,Y)
971 @param size the 2D plan is split into square zones
972 @param length sub path
973 @param max_iter max number of iterations
974 @param pygame pygame for simulation
975 @param screen with pygame
976 @param fLOG logging function
977 @param images save intermediate images
978 @param distance distance function
979 @return path
981 The distance is a function which takes two tuples and returns a distance::
983 def distance(p1, p2):
984 # ...
985 return d
987 Les points identiques sont enlevés puis ajoutés à la fin.
988 """
989 # verification
990 if distance is None:
991 distance = distance_euclidian
992 unique = set()
993 for point in points:
994 if isinstance(point, list):
995 raise TypeError("points cannot be list")
996 unique.add(point)
998 # remove duplicates
999 groups = {}
1000 for p in points:
1001 x, y = p[:2]
1002 if (x, y) in groups:
1003 groups[x, y].append(p)
1004 else:
1005 groups[x, y] = [p]
1007 before = len(points)
1008 points = [v[0] for v in groups.values()]
1009 fLOG("[tsp_kruskal_algorithm] with no duplicates {0} <= {1}".format(
1010 len(points), before))
1012 # begin of the algortihm
1013 fLOG("[tsp_kruskal_algorithm] arbre_poids_minimal nb={0}".format(
1014 len(points)))
1015 di = arbre_poids_minimal(points, size, distance=distance)
1016 arbre = di["arbre"]
1017 X, Y = di["X"], di["Y"]
1018 if screen is not None:
1019 display_arbre(points, arbre, screen=screen, pygame=pygame)
1020 pygame.display.flip()
1021 if images is not None:
1022 c = screen.copy()
1023 for i in range(0, 5):
1024 images.append(c)
1025 fLOG("[tsp_kruskal_algorithm] circuit_eulerien X={0} Y={1}".format(X, Y))
1026 chemin = circuit_eulerien(points, arbre, screen, pygame, fLOG)
1028 if len(chemin) != len(points):
1029 raise Exception("The path should include all points: path:{0} points:{1}".format(
1030 len(chemin), len(points)))
1032 if screen is not None:
1033 display_chemin([points[c] for c in chemin], 0,
1034 screen=screen, pygame=pygame)
1035 pygame.display.flip()
1036 if images is not None:
1037 c = screen.copy()
1038 for i in range(0, 5):
1039 images.append(c)
1041 fLOG("[tsp_kruskal_algorithm] circuit_hamiltonien")
1042 neurone = circuit_hamiltonien(chemin)
1043 neurones = [points[i] for i in neurone]
1044 if screen is not None:
1045 display_chemin(neurones, 0, screen=screen, pygame=pygame)
1046 fLOG("[tsp_kruskal_algorithm] amelioration_chemin")
1047 amelioration_chemin(neurones, size, X, Y, length, screen,
1048 fLOG=fLOG, pygame=pygame, max_iter=max_iter,
1049 images=images, distance=distance)
1051 # we add duplicates back
1052 final = []
1053 for p in neurones:
1054 x, y = p[:2]
1055 g = groups[x, y]
1056 if len(g) == 1:
1057 final.append(p)
1058 else:
1059 final.extend(g)
1060 return final
1063def display_ville(villes, screen, bv, pygame):
1064 """
1065 dessine les villes à l'écran
1066 """
1067 color = 255, 0, 0
1068 color2 = 0, 255, 0
1069 for v in villes:
1070 pygame.draw.circle(screen, color, (int(v[0]), int(v[1])), 3)
1071 v = villes[bv]
1072 pygame.draw.circle(screen, color2, (int(v[0]), int(v[1])), 3)
1075def display_chemin(neurones, bn, screen, pygame):
1076 """
1077 dessine le chemin à l'écran
1078 """
1079 color = 0, 0, 255
1080 color2 = 0, 255, 0
1081 for n in neurones:
1082 pygame.draw.circle(screen, color, (int(n[0]), int(n[1])), 3)
1083 v = neurones[bn]
1084 pygame.draw.circle(screen, color2, (int(v[0]), int(v[1])), 3)
1085 if len(neurones) > 1:
1086 pygame.draw.lines(screen, color, True, neurones, 2)
1089def display_arbre(villes, arbre, mult=1, screen=None, pygame=None):
1090 """
1091 dessine le graphe de poids minimal dꧩni par arbre
1092 """
1093 if mult == 2:
1094 color = 0, 255, 0
1095 li = 4
1096 else:
1097 li = 1
1098 color = 0, 0, 255
1100 for i in range(0, len(villes)):
1101 for j in arbre[i]:
1102 v = (villes[i][0] * mult, villes[i][1] * mult)
1103 vv = (villes[j][0] * mult, villes[j][1] * mult)
1104 pygame.draw.line(screen, color, v, vv, li)
1107def pygame_simulation(size=(800, 500), zone=20, length=10, max_iter=None,
1108 nb=700, fLOG=noLOG, pygame=None, folder=None,
1109 first_click=False, distance=None, flags=0):
1110 """
1111 @param pygame module pygame
1112 @param nb number of cities
1113 @param first_click attend la pression d'un clic de souris avant de commencer
1114 @param folder répertoire où stocker les images de la simulation
1115 @param size taille de l'écran
1116 @param delay delay between two tries
1117 @param folder folder where to save images
1118 @param first_click pause
1119 @param fLOG logging function
1120 @param distance distance function
1121 @param flags see `pygame.display.set_mode <https://www.pygame.org/docs/ref/display.html#pygame.display.set_mode>`_
1122 @return @see fn tsp_kruskal_algorithm
1124 La simulation ressemble à ceci :
1126 .. raw:: html
1128 <video autoplay="" controls="" loop="" height="250">
1129 <source src="http://www.xavierdupre.fr/enseignement/complements/tsp_kruskal.mp4" type="video/mp4" />
1130 </video>
1132 Pour lancer la simulation::
1134 from ensae_teaching_cs.special.tsp_kruskal import pygame_simulation
1135 import pygame
1136 pygame_simulation(pygame)
1138 Voir :ref:`l-tsp_kruskal`.
1139 """
1140 pygame.init()
1141 size = x, y = size[0], size[1]
1142 white = 255, 255, 255
1143 screen = pygame.display.set_mode(size, flags)
1144 screen.fill(white)
1145 points = construit_ville(nb, x, y)
1147 if first_click:
1148 wait_event(pygame)
1150 images = [] if folder is not None else None
1151 tsp_kruskal_algorithm(points, size=zone, length=length, max_iter=max_iter,
1152 fLOG=fLOG, pygame=pygame, screen=screen, images=images,
1153 distance=distance)
1154 fLOG("folder", folder)
1155 fLOG("images", len(images))
1157 if first_click:
1158 wait_event(pygame)
1160 if folder is not None:
1161 fLOG("saving images")
1162 for it, screen in enumerate(images):
1163 if it % 10 == 0:
1164 fLOG("saving image:", it, "/", len(images))
1165 image = os.path.join(folder, "image_%04d.png" % it)
1166 pygame.image.save(screen, image)