Coverage for src/ensae_teaching_cs/special/tsp_kruskal.py: 93%
690 statements
« prev ^ index » next coverage.py v7.1.0, created at 2023-04-28 06:23 +0200
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(f"p{i}: {p[0]},{p[1]}")
349 for i, c in enumerate(chemin):
350 rows.append(f"c{i}: i={c} -> {villes[c][0]},{villes[c][1]}")
351 rows.append(f"bm={bm} ma={ma} bvec={vec2} vec={vec} bl={bl}")
352 rows.append(f"arbre[{bm}]={arbre[bm]}")
353 rows.append(f"arbre[{arbre[bm][0]}]={arbre[arbre[bm][0]]}")
354 mes = "\n".join(rows)
355 raise RuntimeError("this case should not happen\n" + mes)
357 if len(done) < len(villes):
358 # something is wrong (it might an issue with duplicated points)
359 rows = []
360 for i, p in enumerate(villes):
361 rows.append(f"p{i}: {p[0]},{p[1]}")
362 for i, c in enumerate(chemin):
363 rows.append(f"c{i}: i={c} -> {villes[c][0]},{villes[c][1]}")
364 rows.append(f"bm={bm} ma={ma} bvec={vec2} vec={vec} bl={bl}")
365 mes = "\n".join(rows)
366 raise RuntimeError("circuit_eulerien cannot give a path:\n" + mes)
368 return chemin
371def circuit_hamiltonien(chemin):
372 """
373 Extrait un circuit hamiltonien depuis un circuit eulérien,
374 passe par tous les sommets une et une seule fois.
375 """
376 nb = max(chemin) + 1
377 res = []
378 coche = [False for i in range(0, nb)]
379 for c in chemin:
380 if coche[c]:
381 continue
382 res.append(c)
383 coche[c] = True
384 return res
387def equation_droite(p1, p2):
388 """
389 retourne l'équation d'une droite passant par p1 et p2,
390 ax + by + c = 0, retourne les coefficients a,b,c
391 """
392 vec = vecteur_points(p1, p2)
393 a = vec[1]
394 b = - vec[0]
395 c = - a * p1[0] - b * p1[1]
396 return a, b, c
399def evaluation_droite(a, b, c, p):
400 """
401 L'équation d'une droite est : ``ax + by + c``, retourne la valeur
402 de cette expression au point *p*.
403 """
404 return a * p[0] + b * p[1] + c
407def intersection_segment(p1, p2, p3, p4):
408 """
409 Dit si les segments *[p1 p2]* et *[p3 p4]* ont une intersection,
410 ne retourne pas l'intersection.
411 """
412 # équation de la droite (p1 p2)
413 a1, b1, c1 = equation_droite(p1, p2)
414 a2, b2, c2 = equation_droite(p3, p4)
415 s1 = evaluation_droite(a2, b2, c2, p1)
416 s2 = evaluation_droite(a2, b2, c2, p2)
417 s3 = evaluation_droite(a1, b1, c1, p3)
418 s4 = evaluation_droite(a1, b1, c1, p4)
419 return s1 * s2 <= 0 and s3 * s4 <= 0
422def longueur_chemin(chemin, distance):
423 """
424 Retourne la longueur d'un chemin.
425 """
426 s = 0
427 nb = len(chemin)
428 for i in range(0, nb):
429 s += distance(chemin[i], chemin[(i + 1) % nb])
430 return s
433def retournement_essai(chemin, i, j, negligeable=1e-5, distance=None):
434 """
435 Dit s'il est judicieux de parcourir le chemin entre les sommets *i* et *j*
436 en sens inverse, si c'est judicieux, change le chemin et retourne True,
437 sinon, retourne False,
438 si *j < i*, alors la partie à inverser est
439 *i*, *i+1*, *i+2*, ..., *len(chemin)*,
440 *-1*, *0*, *1*, ..., *j*.
441 """
443 nb = len(chemin)
444 if i == j:
445 return False
446 if j - i == -1:
447 return False
448 if j - i - nb == -1:
449 return False
451 ia = (i - 1 + nb) % nb
452 ja = (j + 1) % nb
453 # arcs enlevés
454 d_ia_i = distance(chemin[ia], chemin[i])
455 d_j_ja = distance(chemin[j], chemin[ja])
456 # arcs ajoutés
457 d_ia_j = distance(chemin[ia], chemin[j])
458 d_i_ja = distance(chemin[i], chemin[ja])
459 # amélioration ?
460 d = d_ia_j + d_i_ja - d_j_ja - d_ia_i
461 if d >= -negligeable:
462 return False
464 # si amélioration, il faut retourner le chemin entre les indices i et j
465 jp = j
466 if jp < i:
467 jp = j + nb
468 ip = i
470 while ip < jp:
471 i = ip % nb
472 j = jp % nb
473 ech = chemin[i]
474 chemin[i] = chemin[j]
475 chemin[j] = ech
476 ip = ip + 1
477 jp = jp - 1
479 return True
482def retournement(chemin, taille, fLOG, distance):
483 """
484 Amélioration du chemin par un algorithme simple,
485 utilise des retournements de taille au plus <taille>,
486 retourne le nombre de modifications.
487 """
489 # traitement des petits retournements
490 nb = len(chemin)
491 nb_change = 1
492 nbtout = 0
493 retour = {}
494 while nb_change > 0:
495 nb_change = 0
496 for t in range(1, taille + 1):
497 retour[t] = 0
498 for i in range(0, nb):
499 j = (i + t) % nb
500 b = retournement_essai(chemin, i, j, distance=distance)
501 if b:
502 retour[t] += 1
503 nb_change += 1
504 nbtout += nb_change
505 fLOG("nombre de retournements %d longueur : \t %10.0f --- \t"
506 % (nbtout, longueur_chemin(chemin, distance)), " --- : ", retour)
507 return nbtout
510def echange_position_essai(chemin, a, b, x, inversion, negligeable=1e-5, distance=None):
511 """
512 Echange la place des villes ka et kb pour les placer entre les villes *i* et *i+1*,
513 si inversion est True, on inverse également le chemin inséré, si inversion est False,
514 on ne l'inverse pas,
515 si cela améliore, déplace les villes et retourne True, sinon, retourne False.
516 """
518 nb = len(chemin)
519 xa = (x + 1) % nb
520 ka = (a - 1 + nb) % nb
521 kb = (b + 1) % nb
523 if not inversion:
525 if x == ka:
526 return False
527 if x == kb:
528 return False
529 if xa == ka:
530 return False
531 if b < a:
532 if a <= x <= b + nb:
533 return False
534 elif a <= x <= b:
535 return False
536 if b < a:
537 if a <= x + nb <= b + nb:
538 return False
539 elif a <= x + nb <= b:
540 return False
542 # arcs enlevés
543 d_x_xa = distance(chemin[x], chemin[xa])
544 d_ka_a = distance(chemin[ka], chemin[a])
545 d_b_kb = distance(chemin[b], chemin[kb])
546 # arcs ajoutés
547 d_ka_kb = distance(chemin[ka], chemin[kb])
548 d_x_a = distance(chemin[x], chemin[a])
549 d_b_xa = distance(chemin[b], chemin[xa])
551 d = d_ka_kb + d_x_a + d_b_xa - d_x_xa - d_ka_a - d_b_kb
552 if d >= -negligeable:
553 return False
555 # villes à déplacer
556 ech = []
557 bp = b
558 if bp < a:
559 bp = b + nb
560 for i in range(a, bp + 1):
561 ech.append(chemin[i % nb])
562 diff = bp - a + 1
564 xp = x
565 if xp < b:
566 xp += nb
568 for le in range(b + 1, xp + 1):
569 ll = le % nb
570 bp = (a + le - b - 1) % nb
571 chemin[bp] = chemin[ll]
573 for le in range(0, len(ech)):
574 chemin[(x + le - diff + 1 + nb) % nb] = ech[le]
576 return True
578 else:
580 if x == ka:
581 return False
582 if x == kb:
583 return False
584 if xa == ka:
585 return False
586 if b < a:
587 if a <= x <= b + nb:
588 return False
589 elif a <= x <= b:
590 return False
591 if b < a:
592 if a <= x + nb <= b + nb:
593 return False
594 elif a <= x + nb <= b:
595 return False
597 # arcs enlevés
598 d_x_xa = distance(chemin[x], chemin[xa])
599 d_ka_a = distance(chemin[ka], chemin[a])
600 d_b_kb = distance(chemin[b], chemin[kb])
601 # arcs ajoutés
602 d_ka_kb = distance(chemin[ka], chemin[kb])
603 d_x_b = distance(chemin[x], chemin[b])
604 d_a_xa = distance(chemin[a], chemin[xa])
606 d = d_ka_kb + d_x_b + d_a_xa - d_x_xa - d_ka_a - d_b_kb
607 if d >= -negligeable:
608 return False
610 # villes à déplacer
611 ech = []
612 bp = b
613 if bp < a:
614 bp = b + nb
615 for i in range(a, bp + 1):
616 ech.append(chemin[i % nb])
617 ech.reverse()
618 diff = bp - a + 1
620 xp = x
621 if xp < b:
622 xp += nb
624 for le in range(b + 1, xp + 1):
625 ll = le % nb
626 bp = (a + le - b - 1) % nb
627 chemin[bp] = chemin[ll]
629 for le in range(0, len(ech)):
630 chemin[(x + le - diff + 1 + nb) % nb] = ech[le]
632 return True
635def dessin_arete_zone(chemin, taille_zone, X, Y):
636 """
637 Retourne une liste de listes de listes,
638 ``res[i][j]`` est une liste des arêtes passant près de la zone ``(x,y) = [i][j]``,
639 si *k* in ``res[i][j]``, alors l'arête *k*, *k+1* est dans la zone *(i,j)*,
640 *X* est le nombre de zones horizontalement, *Y* est le nombre de zones verticalement,
641 *taille_zone* est la longueur du côté du carré d'une zone.
642 """
643 res = [[[] for j in range(0, Y + 1)] for i in range(0, X + 1)]
644 nb = len(chemin)
645 mx = min(_[0] for _ in chemin)
646 my = min(_[1] for _ in chemin)
647 for i in range(0, nb):
648 a = chemin[i]
649 b = chemin[(i + 1) % nb]
650 x1, x2 = int(
651 (a[0] - mx) // taille_zone), int((b[0] - mx) // taille_zone)
652 y1, y2 = int(
653 (a[1] - my) // taille_zone), int((b[1] - my) // taille_zone)
654 line = draw_line(x1, y1, x2, y2)
655 for x, y in line:
656 res[x][y].append(i)
657 return res
660def voisinage_zone_xy(x, y, X, Y):
661 """
662 Retourne la liste des voisins d'une zone *(x,y)*
663 sachant qu'il y a *X* zones sur l'axe des abscisses
664 et *Y* zones sur l'axe des ordonnées,
665 inclus *z* dans cette liste
666 """
667 voisin = [(x, y)]
668 if x > 0:
669 voisin.append((x - 1, y))
670 if x < X - 1:
671 voisin.append((x + 1, y))
672 if y > 0:
673 voisin.append((x, y - 1))
674 if y < Y - 1:
675 voisin.append((x, y + 1))
676 if x > 0 and y > 0:
677 voisin.append((x - 1, y - 1))
678 if x > 0 and y < Y - 1:
679 voisin.append((x - 1, y + 1))
680 if x < X - 1 and y > 0:
681 voisin.append((x + 1, y - 1))
682 if x < X - 1 and y < Y - 1:
683 voisin.append((x + 1, y + 1))
684 return voisin
687def echange_position(chemin, taille, taille_zone, X, Y, grande=0.5, fLOG=None, distance=None):
688 """
689 Regarde si on ne peut pas déplacer un segment de longueur taille
690 pour supprimer les arêtes les plus longues,
691 au maximum <grande> longues arêtes,
692 retourne le nombre de changement effectués,
693 *X* est le nombre de zones horizontalement,
694 *Y* est le nombre de zones verticalement,
695 *taille_zone* est la longueur d'un côté du carré d'une zone.
696 """
698 nb = len(chemin)
700 def tri_arete(x, y):
701 """pour trier la liste l par ordre décroissant"""
702 if x[2] < y[2]:
703 return 1
704 elif x[2] > y[2]:
705 return -1
706 else:
707 return 0
709 tmx = min(v[0] for v in chemin)
710 tmy = min(v[1] for v in chemin)
712 # list des arêtes triés par ordre décroissant
713 la = []
714 for i in range(0, nb):
715 im = (i + 1) % nb
716 la.append((i, im, distance(chemin[i], chemin[im])))
717 la = list(sorted(la, key=functools.cmp_to_key(tri_arete)))
719 # zone associée à chaque arête
720 zone = dessin_arete_zone(chemin, taille_zone, X, Y)
722 dseuil = la[int(nb * grande)][2]
723 nbtout = 0
724 nb_change = 0
725 iarete = 0
726 retour = {}
727 for t in range(1, taille + 1):
728 retour[t] = 0
730 while iarete < nb:
731 nb_change = 0
732 arete = la[iarete]
733 iarete += 1
734 x = arete[0]
735 xm = arete[1]
736 a = chemin[x]
737 b = chemin[xm]
738 d = distance(a, b)
739 if d < dseuil:
740 break # arête trop petite
742 # zone traversée par la ligne
743 x1, x2 = (int((a[0] - tmx) // taille_zone),
744 int((b[0] - tmx) // taille_zone))
745 y1, y2 = (int((a[1] - tmy) // taille_zone),
746 int((b[1] - tmy) // taille_zone))
747 ens = draw_line(x1, y1, x2, y2)
748 ville = []
749 for k, l in ens:
750 voisin = voisinage_zone_xy(k, l, X, Y)
751 for u, v in voisin:
752 ville.extend(zone[u][v])
754 # on supprime les doubles
755 ville.sort()
756 if len(ville) == 0:
757 continue
758 sup = []
759 mx = -1
760 for v in ville:
761 if v == mx:
762 sup.append(v)
763 mx = v
764 for s in sup:
765 ville.remove(s)
767 # on étudie les possibilités de casser l'arête (x,xm) aux alentours des villes
768 # comprise dans l'ensemble ville
769 for t in range(1, taille + 1):
771 for i in ville:
773 # on essaye d'insérer le sous-chemin (x- t + 1 + nb) --> x
774 # au milieu de l'arête i,i+1
775 b = echange_position_essai(
776 chemin, (x - t + 1 + nb) % nb, x, i, False, distance=distance)
777 if b:
778 nb_change += 1
779 retour[t] += 1
780 continue
782 # on essaye d'insérer le sous-chemin (xm+ t - 1) --> xm
783 # au milieu de l'arête i,i+1
784 b = echange_position_essai(
785 chemin, (xm + t - 1) % nb, xm, i, False, distance=distance)
786 if b:
787 nb_change += 1
788 retour[t] += 1
789 continue
791 # on essaye de casser l'arête x,xm en insérant
792 # le sous-chemin i --> (i+t) % nb
793 b = echange_position_essai(
794 chemin, i, (i + t) % nb, x, False, distance=distance)
795 if b:
796 nb_change += 1
797 retour[t] += 1
798 continue
799 # idem
800 b = echange_position_essai(
801 chemin, i, (i + t) % nb, x, True, distance=distance)
802 if b:
803 retour[t] += 1
804 nb_change += 1
805 continue
806 # idem
807 b = echange_position_essai(
808 chemin, (i - t + nb) % nb, i, x, False, distance=distance)
809 if b:
810 nb_change += 1
811 retour[t] += 1
812 continue
813 # idem
814 b = echange_position_essai(
815 chemin, (i - t + nb) % nb, i, x, True, distance=distance)
816 if b:
817 retour[t] += 1
818 nb_change += 1
819 continue
821 nbtout += nb_change
823 fLOG("nombre de déplacements %d longueur : \t %10.0f --- \t"
824 % (nbtout, longueur_chemin(chemin, distance=distance)), " --- : ", retour)
825 return nbtout
828def supprime_croisement(chemin, taille_zone, X, Y, fLOG, distance=None):
829 """
830 Supprime les croisements d'arêtes,
831 retourne le nombre de changement effectués,
832 *X* est le nombre de zones horizontalement,
833 *Y* est le nombre de zones verticalement,
834 *taille_zone* est la longueur d'un côté du carré d'une zone
835 """
837 nb = len(chemin)
838 tmx = min(v[0] for v in chemin)
839 tmy = min(v[1] for v in chemin)
841 # zone associée à chaque arête
842 zone = dessin_arete_zone(chemin, taille_zone, X, Y)
843 nbtout = 0
845 for i in range(0, nb):
846 im = (i + 1) % nb
847 a = chemin[i]
848 b = chemin[im]
850 # zone traversée par la ligne
851 x1, x2 = (int((a[0] - tmx) // taille_zone),
852 int((b[0] - tmx) // taille_zone))
853 y1, y2 = (int((a[1] - tmy) // taille_zone),
854 int((b[1] - tmy) // taille_zone))
855 ens = draw_line(x1, y1, x2, y2)
856 ville = []
857 for k, l in ens:
858 voisin = voisinage_zone_xy(k, l, X, Y)
859 for u, v in voisin:
860 ville.extend(zone[u][v])
862 # on supprime les doubles
863 ville.sort()
864 if len(ville) == 0:
865 continue
866 sup = []
867 mx = -1
868 for v in ville:
869 if v == mx:
870 sup.append(v)
871 mx = v
872 for s in sup:
873 ville.remove(s)
875 nb_change = 0
876 for v in ville:
877 b = retournement_essai(chemin, i, v, distance=distance)
878 if b:
879 nb_change += 1
880 continue
881 b = retournement_essai(chemin, im, v, distance=distance)
882 if b:
883 nb_change += 1
884 continue
886 nbtout += nb_change
888 fLOG("nombre de croisements %d longueur : \t %10.0f --- \t"
889 % (nbtout, longueur_chemin(chemin, distance=distance)))
890 return nbtout
893def amelioration_chemin(chemin, taille_zone, X, Y, taille=10, screen=None,
894 fLOG=None, pygame=None, max_iter=None, images=None,
895 distance=None):
896 """
897 Amélioration du chemin par un algorithme simple,
898 utilise des retournements de taille au plus *taille*,
899 traite les arcs qui se croisent,
900 traite les grands arcs, utilise un quadrillage de taille *window*,
901 *X* est le nombre de zones horizontalement,
902 *Y* est le nombre de zones verticalement,
903 *taille_zone* est la longueur d'un côté du carré d'une zone.
904 """
906 white = 255, 255, 255
908 if pygame is not None and images is not None:
909 images.append(screen.copy())
911 # première étape rapide
912 iter = 0
913 nb = 1
914 while nb > 0 and (max_iter is None or iter < max_iter):
915 nb = retournement(chemin, taille, fLOG=fLOG, distance=distance)
916 if screen is not None:
917 screen.fill(white)
918 display_chemin(chemin, 0, screen, pygame=pygame)
919 pygame.display.flip()
920 if images is not None:
921 images.append(screen.copy())
922 empty_main_loop(pygame)
923 iter += 1
925 # amélioration
926 nb = 1
927 while nb > 0 and (max_iter is None or iter < max_iter):
928 nb = retournement(chemin, taille, fLOG=fLOG, distance=distance)
929 if screen is not None:
930 screen.fill(white)
931 display_chemin(chemin, 0, screen=screen, pygame=pygame)
932 pygame.display.flip()
933 if images is not None:
934 images.append(screen.copy())
935 empty_main_loop(pygame)
936 nb += echange_position(chemin, taille // 2,
937 taille_zone, X, Y, fLOG=fLOG,
938 distance=distance)
939 if screen is not None:
940 screen.fill(white)
941 display_chemin(chemin, 0, screen=screen, pygame=pygame)
942 pygame.display.flip()
943 if images is not None:
944 images.append(screen.copy())
945 empty_main_loop(pygame)
946 nb += supprime_croisement(chemin, taille_zone,
947 X, Y, fLOG=fLOG, distance=distance)
948 if screen is not None:
949 screen.fill(white)
950 display_chemin(chemin, 0, screen=screen, pygame=pygame)
951 pygame.display.flip()
952 if images is not None:
953 images.append(screen.copy())
954 empty_main_loop(pygame)
955 iter += 1
958def tsp_kruskal_algorithm(points, size=20, length=10, max_iter=None,
959 fLOG=noLOG, pygame=None, screen=None, images=None,
960 distance=None):
961 """
962 Finds the shortest path going through all points,
963 points require to be a 2 dimensional space.
965 @param points list 2-tuple (X,Y)
966 @param size the 2D plan is split into square zones
967 @param length sub path
968 @param max_iter max number of iterations
969 @param pygame pygame for simulation
970 @param screen with pygame
971 @param fLOG logging function
972 @param images save intermediate images
973 @param distance distance function
974 @return path
976 The distance is a function which takes two tuples and returns a distance::
978 def distance(p1, p2):
979 # ...
980 return d
982 Les points identiques sont enlevés puis ajoutés à la fin.
983 """
984 # verification
985 if distance is None:
986 distance = distance_euclidian
987 unique = set()
988 for point in points:
989 if isinstance(point, list):
990 raise TypeError("points cannot be list")
991 unique.add(point)
993 # remove duplicates
994 groups = {}
995 for p in points:
996 x, y = p[:2]
997 if (x, y) in groups:
998 groups[x, y].append(p)
999 else:
1000 groups[x, y] = [p]
1002 before = len(points)
1003 points = [v[0] for v in groups.values()]
1004 fLOG(
1005 f"[tsp_kruskal_algorithm] with no duplicates {len(points)} <= {before}")
1007 # begin of the algortihm
1008 fLOG(f"[tsp_kruskal_algorithm] arbre_poids_minimal nb={len(points)}")
1009 di = arbre_poids_minimal(points, size, distance=distance)
1010 arbre = di["arbre"]
1011 X, Y = di["X"], di["Y"]
1012 if screen is not None:
1013 display_arbre(points, arbre, screen=screen, pygame=pygame)
1014 pygame.display.flip()
1015 if images is not None:
1016 c = screen.copy()
1017 for i in range(0, 5):
1018 images.append(c)
1019 fLOG(f"[tsp_kruskal_algorithm] circuit_eulerien X={X} Y={Y}")
1020 chemin = circuit_eulerien(points, arbre, screen, pygame, fLOG)
1022 if len(chemin) != len(points):
1023 raise RuntimeError( # pragma: no cover
1024 "The path should include all points: path:{0} points:{1}".format(
1025 len(chemin), len(points)))
1027 if screen is not None:
1028 display_chemin([points[c] for c in chemin], 0,
1029 screen=screen, pygame=pygame)
1030 pygame.display.flip()
1031 if images is not None:
1032 c = screen.copy()
1033 for i in range(0, 5):
1034 images.append(c)
1036 fLOG("[tsp_kruskal_algorithm] circuit_hamiltonien")
1037 neurone = circuit_hamiltonien(chemin)
1038 neurones = [points[i] for i in neurone]
1039 if screen is not None:
1040 display_chemin(neurones, 0, screen=screen, pygame=pygame)
1041 fLOG("[tsp_kruskal_algorithm] amelioration_chemin")
1042 amelioration_chemin(neurones, size, X, Y, length, screen,
1043 fLOG=fLOG, pygame=pygame, max_iter=max_iter,
1044 images=images, distance=distance)
1046 # we add duplicates back
1047 final = []
1048 for p in neurones:
1049 x, y = p[:2]
1050 g = groups[x, y]
1051 if len(g) == 1:
1052 final.append(p)
1053 else:
1054 final.extend(g)
1055 return final
1058def display_ville(villes, screen, bv, pygame):
1059 """
1060 dessine les villes à l'écran
1061 """
1062 color = 255, 0, 0
1063 color2 = 0, 255, 0
1064 for v in villes:
1065 pygame.draw.circle(screen, color, (int(v[0]), int(v[1])), 3)
1066 v = villes[bv]
1067 pygame.draw.circle(screen, color2, (int(v[0]), int(v[1])), 3)
1070def display_chemin(neurones, bn, screen, pygame):
1071 """
1072 dessine le chemin à l'écran
1073 """
1074 color = 0, 0, 255
1075 color2 = 0, 255, 0
1076 for n in neurones:
1077 pygame.draw.circle(screen, color, (int(n[0]), int(n[1])), 3)
1078 v = neurones[bn]
1079 pygame.draw.circle(screen, color2, (int(v[0]), int(v[1])), 3)
1080 if len(neurones) > 1:
1081 pygame.draw.lines(screen, color, True, neurones, 2)
1084def display_arbre(villes, arbre, mult=1, screen=None, pygame=None):
1085 """
1086 dessine le graphe de poids minimal dꧩni par arbre
1087 """
1088 if mult == 2:
1089 color = 0, 255, 0
1090 li = 4
1091 else:
1092 li = 1
1093 color = 0, 0, 255
1095 for i in range(0, len(villes)):
1096 for j in arbre[i]:
1097 v = (villes[i][0] * mult, villes[i][1] * mult)
1098 vv = (villes[j][0] * mult, villes[j][1] * mult)
1099 pygame.draw.line(screen, color, v, vv, li)
1102def pygame_simulation(size=(800, 500), zone=20, length=10, max_iter=None,
1103 nb=700, fLOG=noLOG, pygame=None, folder=None,
1104 first_click=False, distance=None, flags=0):
1105 """
1106 @param pygame module pygame
1107 @param nb number of cities
1108 @param first_click attend la pression d'un clic de souris avant de commencer
1109 @param folder répertoire où stocker les images de la simulation
1110 @param size taille de l'écran
1111 @param delay delay between two tries
1112 @param folder folder where to save images
1113 @param first_click pause
1114 @param fLOG logging function
1115 @param distance distance function
1116 @param flags see `pygame.display.set_mode <https://www.pygame.org/docs/ref/display.html#pygame.display.set_mode>`_
1117 @return @see fn tsp_kruskal_algorithm
1119 La simulation ressemble à ceci :
1121 .. raw:: html
1123 <video autoplay="" controls="" loop="" height="250">
1124 <source src="http://www.xavierdupre.fr/enseignement/complements/tsp_kruskal.mp4" type="video/mp4" />
1125 </video>
1127 Pour lancer la simulation::
1129 from ensae_teaching_cs.special.tsp_kruskal import pygame_simulation
1130 import pygame
1131 pygame_simulation(pygame)
1133 Voir :ref:`l-tsp_kruskal`.
1134 """
1135 pygame.init()
1136 size = x, y = size[0], size[1]
1137 white = 255, 255, 255
1138 screen = pygame.display.set_mode(size, flags)
1139 screen.fill(white)
1140 points = construit_ville(nb, x, y)
1142 if first_click:
1143 wait_event(pygame) # pragma: no cover
1145 images = [] if folder is not None else None
1146 tsp_kruskal_algorithm(points, size=zone, length=length, max_iter=max_iter,
1147 fLOG=fLOG, pygame=pygame, screen=screen, images=images,
1148 distance=distance)
1149 fLOG("folder", folder)
1150 fLOG("images", len(images))
1152 if first_click:
1153 wait_event(pygame) # pragma: no cover
1155 if folder is not None:
1156 fLOG("saving images")
1157 for it, screen in enumerate(images):
1158 if it % 10 == 0:
1159 fLOG("saving image:", it, "/", len(images))
1160 image = os.path.join(folder, "image_%04d.png" % it)
1161 pygame.image.save(screen, image)