Coverage for src/ensae_projects/challenge/city_tour.py : 86%
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"""
2@file
3@brief Function to solve the problem of the
4`Route Inspection Problem <https://en.wikipedia.org/wiki/Route_inspection_problem>`_.
5"""
6import math
7import random
8import fractions
11class SolutionException(Exception):
12 """
13 wrong solution
14 """
15 pass
18def haversine_distance(lat1, lng1, lat2, lng2):
19 """
20 Computes `Haversine formula <http://en.wikipedia.org/wiki/Haversine_formula>`_.
22 .. index:: Haversine
24 @param lat1 lattitude
25 @param lng1 longitude
26 @param lat2 lattitude
27 @param lng2 longitude
28 @return distance
29 """
30 radius = 6371
31 dlat = math.radians(lat2 - lat1)
32 dlon = math.radians(lng2 - lng1)
33 a = math.sin(dlat / 2) * math.sin(dlat / 2) + math.cos(math.radians(lat1)) \
34 * math.cos(math.radians(lat2)) * math.sin(dlon / 2) * math.sin(dlon / 2)
35 c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
36 d = radius * c
37 return d
40def distance_solution(edges_index, edges, distances, solution, exc=True):
41 """
42 Checks if a solution is really a solution and returns the distance of it,
43 None if it is not a solution. The function does not case about the order.
45 @param edges_index list of indices of edges (if None --> range(len(edges))
46 @param edges list of tuple (vertex A, vertex B)
47 @param distances list of distances of each edge
48 @param solution proposed solutions (list of edge indices)
49 @param exc raises an exception in case of error
50 """
51 if edges_index is None:
52 edges_index = list(range(len(edges)))
53 indices = set(edges_index)
54 solset = set(solution)
55 if len(indices) != len(solset):
56 if exc:
57 raise SolutionException("Different number of distinct edges:\nexpected={0} got={1}\n"
58 "Did you cover all the edges?".format(len(indices), len(solset)))
59 return None
61 for s in solution:
62 if s not in indices:
63 raise SolutionException(
64 "Index {0} is not in edges_index".format(s))
66 doubles = {}
67 for a, b in edges:
68 if a > b:
69 a, b = b, a
70 if (a, b) in doubles:
71 if exc:
72 raise SolutionException(
73 "Edge {0} is duplicated in edges.".format((a, b)))
74 return None
75 doubles[a, b] = 1
77 corres = {e: i for i, e in enumerate(edges_index)}
78 degrees = {}
79 for s in solution:
80 a, b = edges[corres[s]]
81 degrees[a] = degrees.get(a, 0) + 1
82 degrees[b] = degrees.get(b, 0) + 1
84 odd, even = 0, 0
85 for d in degrees.values():
86 if d % 2 == 0:
87 even += 1
88 else:
89 odd += 1
91 if odd > 2:
92 if exc:
93 red = list(sorted([(k, v)
94 for k, v in degrees.items() if v % 2 != 0]))
95 if len(red) > 10:
96 red = red[:10] + ["..."]
97 raise SolutionException(
98 "Are you sure? The path is inconsistent. Some help:\n" + str(red))
99 return None
100 else:
101 return sum(distances[corres[s]] for s in solution)
104def compute_degrees(edges):
105 """
106 Compute the degree of vertices.
108 @param edges list of tuple
109 @return dictionary {key: degree}
110 """
111 res = {}
112 for a, b in edges:
113 res[a] = res.get(a, 0) + 1
114 res[b] = res.get(b, 0) + 1
115 return res
118def euler_path(edges_index, edges, solution):
119 """
120 Computes an Eulerian path.
122 .. index:: Eulerian
124 @param edges_index list of indices of edges (if None --> range(len(edges))
125 @param edges list of tuple (vertex A, vertex B)
126 @param solution proposed solutions (list of edge indices)
127 @return path, list of edges indices
129 The function assumes every vertex in the graph defined by *edges*
130 has an even degree.
131 """
132 if edges_index is None:
133 edges_index = range(len(edges))
134 pos = {k: i for i, k in enumerate(edges_index)}
135 indices = [pos[s] for s in solution]
136 edges = [edges[i] for i in indices]
138 degrees = compute_degrees(edges)
139 odd = {k: v for k, v in degrees.items() if v % 2 == 1}
140 if len(odd) not in (0, 2):
141 odd = list(sorted((k, v) for k, v in odd.items()))
142 if len(odd) > 10:
143 odd = odd[:10] + ["..."]
144 raise SolutionException(
145 "Some vertices have an odd degree. This is not allowed.\n" + str(odd))
147 if len(odd) == 2:
148 # we add an extra edge which we remove later
149 odd = list(odd.keys())
150 remove = (odd[0], odd[1])
151 edges.append(remove)
152 else:
153 remove = None
155 res = _euler_path(edges)
156 pathi = [_[1][0] for _ in res]
157 if remove is not None:
158 index = pathi.index(len(edges) - 1)
159 if index == 0:
160 pathi = pathi[1:]
161 elif index == len(edges) - 1:
162 pathi = pathi[:-1]
163 else:
164 pathi = pathi[index + 1:] + pathi[:index]
165 pathi = [solution[i] for i in pathi]
166 return pathi
169def _euler_path(edges):
170 """
171 Computes an Eulerian path.
173 @param edges edges
174 @return path, list of (vertex,edge)
176 The function assumes every vertex in the graph defined by *edges*
177 has an even degree.
178 """
179 alledges = {}
180 edges_from = {}
181 for i, k in enumerate(edges):
182 if isinstance(k, list):
183 k = tuple(k)
184 v = (i,) + k
185 alledges[k] = v
186 a, b = k
187 alledges[b, a] = alledges[a, b]
188 if a not in edges_from:
189 edges_from[a] = []
190 if b not in edges_from:
191 edges_from[b] = []
192 edges_from[a].append(alledges[a, b])
193 edges_from[b].append(alledges[a, b])
195 degre = {}
196 for a, v in edges_from.items():
197 t = len(v)
198 degre[t] = degre.get(t, 0) + 1
200 two = [a for a, v in edges_from.items() if len(v) == 2]
201 odd = [a for a, v in edges_from.items() if len(v) % 2 == 1]
202 if len(odd) not in (0, 2):
203 add = "\n" + str(odd) if len(odd) < 10 else ""
204 raise SolutionException(
205 "Some vertices have an odd degree. This is not allowed." + add)
207 begin = two[0]
209 # checking
210 for v, le in edges_from.items():
211 for e in le:
212 to = e[1] if v != e[1] else e[2]
213 if to not in edges_from:
214 raise SolutionException(
215 "Unable to find vertex {0} for edge ({0},{1}).".format(to, v))
216 if to == v:
217 raise SolutionException("Circular edge {0}.".format(to))
219 # loop
220 path = _explore_path(edges_from, begin)
221 for p in path:
222 if len(p) == 0:
223 raise NotImplementedError("This exception should not happen.")
224 while len(edges_from) > 0:
225 start = None
226 for i, p in enumerate(path):
227 if p[0] in edges_from:
228 start = i, p
229 break
230 sub = _explore_path(edges_from, start[1][0])
231 i = start[0]
232 path[i:i + 1] = path[i:i + 1] + sub
233 return path
236def _delete_edge(edges_from, n, to):
237 """
238 Removes an edge from the graph.
240 @param edges_from structure which contains the edges (will be modified)
241 @param n first vertex
242 @param to second vertex
243 @return the edge
244 """
245 le = edges_from[to]
246 f = None
247 for i, e in enumerate(le):
248 if (e[1] == to and e[2] == n) or (e[2] == to and e[1] == n):
249 f = i
250 break
252 assert f is not None
253 del le[f]
254 if len(le) == 0:
255 del edges_from[to]
257 le = edges_from[n]
258 f = None
259 for i, e in enumerate(le):
260 if (e[1] == to and e[2] == n) or (e[2] == to and e[1] == n):
261 f = i
262 break
264 assert f is not None
265 keep = le[f]
266 del le[f]
267 if len(le) == 0:
268 del edges_from[n]
270 return keep
273def _explore_path(edges_from, begin):
274 """
275 Explores an Eulerian path, remove used edges from *edges_from*.
277 @param edges_from structure which contains the edges (will be modified)
278 @param begin first vertex to use
279 @return path
280 """
281 path = [(begin, None)]
282 stay = True
283 while stay and len(edges_from) > 0:
285 n = path[-1][0]
286 if n not in edges_from:
287 # fin
288 break
289 le = edges_from[n]
291 if len(le) == 1:
292 h = 0
293 e = le[h]
294 to = e[1] if n != e[1] else e[2]
295 else:
296 to = None
297 nb = 100
298 while to is None or to == begin:
299 h = random.randint(0, len(le) - 1) if len(le) > 1 else 0
300 e = le[h]
301 to = e[1] if n != e[1] else e[2]
302 nb -= 1
303 if nb < 0:
304 raise NotImplementedError(
305 "Algorithm issue {0}".format(len(path)))
307 if len(edges_from[to]) == 1:
308 if begin != to:
309 raise NotImplementedError("Wrong algorithm.")
310 stay = False
312 keep = _delete_edge(edges_from, n, to)
313 path.append((to, keep))
315 return path[1:]
318def distance_vertices(edges, vertices, distances):
319 """
320 Computes the length of edges if distances is None.
322 @param edges list of tuple (vertex A, vertex B)
323 @param vertices locations of the vertices
324 @param distances distances (None or list of floats)
325 @return distances (list of float)
326 """
327 if distances is None:
328 distances = []
329 while len(distances) < len(edges):
330 distances.append(None)
331 for i, edge in enumerate(edges):
332 if distances[i] is not None:
333 continue
334 a, b = edge
335 va = vertices[a]
336 vb = vertices[b]
337 d = haversine_distance(va[0], va[1], vb[0], vb[1])
338 distances[i] = d
339 return distances
342def bellman_distances(edges, distances, fLOG=None):
343 """
344 Computes shortest distances between all vertices.
345 It assumes edges are symmetric.
347 @param edges list of tuple (vertex A, vertex B)
348 @param distances distances (list of floats)
349 @param fLOG logging function
350 @return dictionary of distances
352 This function could be implemented based on
353 `shortest_path <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csgraph.shortest_path.html>`_.
354 """
355 dist = {(a, b): d for d, (a, b) in zip(distances, edges)}
356 dist.update({(b, a): d for d, (a, b) in zip(distances, edges)})
357 dist0 = dist.copy()
359 iter = 0 # pylint: disable=W0622
360 up = dist
361 while len(up) > 0:
362 iter += 1
363 up = {}
364 for (a, b), d1 in dist.items():
365 if (a, b) not in dist0 and (a, b) not in up:
366 # no need to continue as no update during the last iteration
367 continue
368 for (aa, bb), d2 in dist.items():
369 # not the most efficient
370 if b == aa and a != bb:
371 d = d1 + d2
372 if (a, bb) not in dist or dist[a, bb] > d:
373 up[a, bb] = d
374 up[bb, a] = d
375 dist.update(up)
376 if fLOG:
377 sum_values = sum(dist.values())
378 mes = "iteration={0} modif={1} #dist={2} sum_values={3} avg={4}".format(
379 iter, len(up), len(dist), sum_values, sum_values / len(dist))
380 fLOG("[bellman_distances] " + mes)
382 return dist
385def dijkstra_path(edges, distances, va, vb):
386 """
387 Returns the best path between two vertices.
388 Uses `Dikjstra <https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm>`_ algorithm.
390 @param edges list of edges.
391 @param distances list of distances
392 @param va first vertex
393 @param vb last vertex
394 @return list of edges
396 This function could be implemented based on
397 `shortest_path <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csgraph.shortest_path.html>`_.
398 """
399 dist = {va: 0}
400 prev = {va: None}
401 modif = 1
402 while modif > 0:
403 modif = 0
404 for (a, b), d in zip(edges, distances):
405 if a in dist:
406 d2 = dist[a] + d
407 if b not in dist or dist[b] > d2:
408 dist[b] = d2
409 prev[b] = a
410 modif += 1
411 if b in dist:
412 d2 = dist[b] + d
413 if a not in dist or dist[a] > d2:
414 dist[a] = d2
415 prev[a] = b
416 modif += 1
417 rev = {(a, b): i for i, (a, b) in enumerate(edges)}
418 rev.update({(b, a): i for i, (a, b) in enumerate(edges)})
419 path = []
420 v = vb
421 while v is not None:
422 path.append(v)
423 v = prev[v]
424 path.reverse()
425 return [rev[a, b] for a, b in zip(path[:-1], path[1:])]
428def matching_vertices(distances, algo="blossom"):
429 """
430 Finds the best match between vertices.
432 @param distances result of function @bellman_distances but
433 only for odd vertices (odd = odd degree),
434 dictionary { (vi,vj) : distance}
435 @param algo algorithm (see below)
436 @return sequences of best matches.
438 If ``algo=='hungarian'``,
439 the function relies on `linear_sum_assignment
440 <https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.linear_sum_assignment.html>`_
441 from `scipy <https://docs.scipy.org/doc/scipy-0.18.1/>`_ which uses the
442 `Hungarian Algorithm <https://en.wikipedia.org/wiki/Hungarian_algorithm>`_.
443 Vertex index do not have to start from zero and be continuous.
444 The function will handle that particular case. However, the method
445 works for a bi-partite matching and is not suitable here unless
446 the algorithm is modified. Not done (yet?).
448 If ``algo=='blossom'``, it uses the
449 `Blossom algorithm <https://en.wikipedia.org/wiki/Blossom_algorithm>`_
450 which is known to be in :math:`O(n^3)` and finds the optimal matching.
452 If ``algo=='basic'``, the function sorts the distances
453 by increasing order and builds new edges as they come in this list.
454 It does not return an optimal solution but is much faster
455 when the graph is big.
456 """
457 if not isinstance(distances, dict):
458 raise TypeError(
459 "Unexpected type for distances, this should be a dictionary.")
460 from numpy import empty
461 from scipy.optimize import linear_sum_assignment
462 unique = set(_[0] for _ in distances) | set(_[1] for _ in distances)
463 mapping = {}
464 rev = {}
465 for s in unique:
466 n = len(mapping)
467 mapping[s] = n
468 rev[n] = s
469 mx = len(mapping)
471 # Hungarian algorithm
472 if algo == "hungarian":
473 cost = empty((mx, mx))
474 mx = abs(max(distances.values())) * len(mapping) * 10
475 cost[:] = mx
476 for (i, j), v in distances.items():
477 if i == j:
478 raise ValueError(
479 "Unreasonable case: {0} == {1}, v={2}".format(i, j, v))
480 cost[mapping[i], mapping[j]] = v
482 row_ind, col_ind = linear_sum_assignment( # pylint: disable=W0632
483 cost) # pylint: disable=W0632
484 pairs = [(rev[i], rev[j]) for i, j in zip(row_ind, col_ind)]
485 for a, b in pairs:
486 if a == b:
487 raise ValueError("Issue with one pair a == b.\n{0}".format(
488 "\n".join(str(_) for _ in pairs)))
490 # we remove duplicates
491 done = set()
492 final = []
493 for a, b in pairs:
494 if (a, b) not in done:
495 final.append((a, b))
496 done.add((b, a))
497 if len(final) * 2 != len(pairs):
498 mes = "final={0}\n{3}={1}\ncost\n{2}".format(
499 final, pairs, cost, algo)
500 raise ValueError(
501 "Did you use the tweak? The matching should be symmetric.\n" + mes)
502 return final
503 elif algo == "basic":
504 # we sort pair by increasing order
505 dists = [(v, i, j) for (i, j), v in distances.items()]
506 dists.sort()
508 # the graph is fully connected so we don't need to check
509 # if selecting an edge leads to an impossible solution
510 vdone = set()
511 pairs = []
512 for v, i, j in dists:
513 if i not in vdone and j not in vdone:
514 pairs.append((i, j))
515 vdone.add(i)
516 vdone.add(j)
517 return pairs
519 elif algo == "blossom":
520 from .blossom import Vertex, StructureUpToDate, TreeStructureChanged, MaximumDualReached, INF
522 vertices = {}
523 for v in mapping:
524 vertices[v] = Vertex(v)
525 for (i, j), v in distances.items():
526 f = fractions.Fraction(v)
527 vertices[i].add_edge_to(vertices[j], f)
529 def update_tree_structures(roots):
530 try:
531 while True:
532 try:
533 for root in roots:
534 root.alter_tree(roots)
535 raise StructureUpToDate()
536 except TreeStructureChanged:
537 pass
538 except StructureUpToDate:
539 pass
541 def get_max_delta(roots):
542 if len(roots) == 0:
543 raise MaximumDualReached()
544 delta = INF
545 for root in roots:
546 delta = min(delta, root.get_max_delta())
547 assert delta >= 0
548 return delta
550 roots = set(vertices.values())
552 try:
553 while True:
554 delta = get_max_delta(roots)
555 for root in roots:
556 root.adjust_charge(delta)
557 update_tree_structures(roots)
558 except MaximumDualReached:
559 # done
560 pass
562 M = set()
563 for v in vertices.values():
564 M.update(e for e in v.edges if e.selected)
566 pairs = [e.extremities for e in M]
567 return pairs
568 else:
569 raise NotImplementedError("Not recognized: {0}".format(algo))
572def best_euler_path(edges, vertices, distances, edges_index=None, algo="blossom", fLOG=None):
573 """
574 Computes the final solution for the Eulerian path.
576 @param edges edges (list of tuple)
577 @param vertices location of the vertices (not needed if distances are filled)
578 @param distances distances for each edge (list of distance)
579 @param edges_index list of indices of edges (if None --> range(len(edges))
580 @param algo algorithm to use to computes the matching
581 (see @see fn matching_vertices)
582 @param fLOG logging function
583 @return list of edges as tuple, list of edges as indices, distance
584 """
585 if fLOG:
586 fLOG("[best_euler_path] distance_vertices #edges={0}".format(
587 len(edges)))
588 distances = distance_vertices(edges, vertices, distances)
589 if fLOG:
590 fLOG("[best_euler_path] bellman_distances")
591 dist = bellman_distances(edges, distances, fLOG=fLOG)
592 if fLOG:
593 fLOG("[best_euler_path] degrees and distances between odd vertices")
594 degrees = compute_degrees(edges)
595 odd = {k: v for k, v in degrees.items() if v % 2 != 0}
596 odd_dist = {k: v for k, v in dist.items() if k[0] in odd and k[1] in odd}
597 if fLOG:
598 fLOG("[best_euler_path] matching #odd={0}, #odd_dist={1}".format(
599 len(odd), len(odd_dist)))
600 pairs = matching_vertices(odd_dist, algo=algo)
601 if fLOG:
602 fLOG("[best_euler_path] build solution")
603 solution = list(edges)
604 for va, vb in pairs:
605 dij = dijkstra_path(edges, distances, va, vb)
606 solution.extend([edges[e] for e in dij])
607 if fLOG:
608 fLOG("[best_euler_path] order edges to get the path #edges={0}".format(
609 len(solution)))
610 if edges_index is None:
611 edges_index = list(range(len(edges)))
612 mapping = {}
613 rev = {}
614 for edge, index in zip(edges, edges_index):
615 mapping[edge] = index
616 rev[index] = edge
617 sol_index = [mapping[e] for e in solution]
618 euler_index = euler_path(edges_index, edges, sol_index)
619 euler = [rev[i] for i in euler_index]
620 if fLOG:
621 fLOG("[best_euler_path] done.")
622 d = distance_solution(edges_index, edges, distances, sol_index)
623 return euler, euler_index, d