Coverage for src/ensae_teaching_cs/special/tsp_kruskal.py: 93%

690 statements  

« prev     ^ index     » next       coverage.py v7.1.0, created at 2023-01-27 05:44 +0100

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 

17 

18 

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 

44 

45 

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) 

53 

54 

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]) 

60 

61 

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]) 

67 

68 

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) 

82 

83 

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) 

97 

98 

99def oppose_vecteur(vec): 

100 """ 

101 retourne le vecteur opposé. 

102 """ 

103 return (-vec[0], -vec[1]) 

104 

105 

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 : 

110 

111 - les coordonnées de la ville 

112 - ses coordonnées en zone, (zx, zy) 

113 - son indice dans la liste villes 

114 

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 

126 

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)) 

137 

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])) 

143 

144 if ask_zone: 

145 return zones, X, Y, mx, my, Zmax, zone 

146 return zones, X, Y, mx, my, Zmax 

147 

148 

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 

178 

179 

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. 

187 

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 """ 

196 

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 

204 

205 rz = repartition_zone(villes, zone_taille=zone_taille) 

206 zones, X, Y, mx, my, Zmax = rz[:6] 

207 

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)) 

217 

218 # tri 

219 li = list(sorted(li, key=functools.cmp_to_key(tri_distance))) 

220 

221 # nombre de composantes connexes 

222 nb_comp = len(villes) 

223 

224 # indice de la composante d'une ville 

225 num_comp = [i for i in range(0, len(villes))] 

226 

227 # liste des voisins pour chaque ville 

228 arbre = [[] for i in range(0, len(villes))] 

229 

230 # liste des villes par composante connexe 

231 list_comp = [[i] for i in range(0, len(villes))] 

232 

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 

251 

252 if nb_comp <= 1: 

253 break # il n'y a plus qu'une seule composante connexe, inutile de continuer 

254 

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])] 

268 

269 return dict(arbre=arbre, X=X, Y=Y, mx=mx, my=my, Zmax=Zmax) 

270 

271 

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 """ 

281 

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) 

289 

290 bm = random.randint(0, len(has) - 1) 

291 bm = has[bm] 

292 

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 

336 

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 Exception("this case should not happen\n" + mes) 

356 

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 Exception("circuit_eulerien cannot give a path:\n" + mes) 

367 

368 return chemin 

369 

370 

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 

385 

386 

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 

397 

398 

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 

405 

406 

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 

420 

421 

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 

431 

432 

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 """ 

442 

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 

450 

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 

463 

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 

469 

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 

478 

479 return True 

480 

481 

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 """ 

488 

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 

508 

509 

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 """ 

517 

518 nb = len(chemin) 

519 xa = (x + 1) % nb 

520 ka = (a - 1 + nb) % nb 

521 kb = (b + 1) % nb 

522 

523 if not inversion: 

524 

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 

541 

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]) 

550 

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 

554 

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 

563 

564 xp = x 

565 if xp < b: 

566 xp += nb 

567 

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] 

572 

573 for le in range(0, len(ech)): 

574 chemin[(x + le - diff + 1 + nb) % nb] = ech[le] 

575 

576 return True 

577 

578 else: 

579 

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 

596 

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]) 

605 

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 

609 

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 

619 

620 xp = x 

621 if xp < b: 

622 xp += nb 

623 

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] 

628 

629 for le in range(0, len(ech)): 

630 chemin[(x + le - diff + 1 + nb) % nb] = ech[le] 

631 

632 return True 

633 

634 

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 

658 

659 

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 

685 

686 

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 """ 

697 

698 nb = len(chemin) 

699 

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 

708 

709 tmx = min(v[0] for v in chemin) 

710 tmy = min(v[1] for v in chemin) 

711 

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))) 

718 

719 # zone associée à chaque arête 

720 zone = dessin_arete_zone(chemin, taille_zone, X, Y) 

721 

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 

729 

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 

741 

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]) 

753 

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) 

766 

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): 

770 

771 for i in ville: 

772 

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 

781 

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 

790 

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 

820 

821 nbtout += nb_change 

822 

823 fLOG("nombre de déplacements %d longueur : \t %10.0f --- \t" 

824 % (nbtout, longueur_chemin(chemin, distance=distance)), " --- : ", retour) 

825 return nbtout 

826 

827 

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 """ 

836 

837 nb = len(chemin) 

838 tmx = min(v[0] for v in chemin) 

839 tmy = min(v[1] for v in chemin) 

840 

841 # zone associée à chaque arête 

842 zone = dessin_arete_zone(chemin, taille_zone, X, Y) 

843 nbtout = 0 

844 

845 for i in range(0, nb): 

846 im = (i + 1) % nb 

847 a = chemin[i] 

848 b = chemin[im] 

849 

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]) 

861 

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) 

874 

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 

885 

886 nbtout += nb_change 

887 

888 fLOG("nombre de croisements %d longueur : \t %10.0f --- \t" 

889 % (nbtout, longueur_chemin(chemin, distance=distance))) 

890 return nbtout 

891 

892 

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 """ 

905 

906 white = 255, 255, 255 

907 

908 if pygame is not None and images is not None: 

909 images.append(screen.copy()) 

910 

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 

924 

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 

956 

957 

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. 

964 

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 

975 

976 The distance is a function which takes two tuples and returns a distance:: 

977 

978 def distance(p1, p2): 

979 # ... 

980 return d 

981 

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) 

992 

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] 

1001 

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}") 

1006 

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) 

1021 

1022 if len(chemin) != len(points): 

1023 raise Exception( # pragma: no cover 

1024 "The path should include all points: path:{0} points:{1}".format( 

1025 len(chemin), len(points))) 

1026 

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) 

1035 

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) 

1045 

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 

1056 

1057 

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) 

1068 

1069 

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) 

1082 

1083 

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 

1094 

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) 

1100 

1101 

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 

1118 

1119 La simulation ressemble à ceci : 

1120 

1121 .. raw:: html 

1122 

1123 <video autoplay="" controls="" loop="" height="250"> 

1124 <source src="http://www.xavierdupre.fr/enseignement/complements/tsp_kruskal.mp4" type="video/mp4" /> 

1125 </video> 

1126 

1127 Pour lancer la simulation:: 

1128 

1129 from ensae_teaching_cs.special.tsp_kruskal import pygame_simulation 

1130 import pygame 

1131 pygame_simulation(pygame) 

1132 

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) 

1141 

1142 if first_click: 

1143 wait_event(pygame) # pragma: no cover 

1144 

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)) 

1151 

1152 if first_click: 

1153 wait_event(pygame) # pragma: no cover 

1154 

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)