# Solving the chinese postman problem

Some pre-Halloween post today. It started actually while I was in Barcelona : kids wanted to go back to some store we’ve seen the first day, in the gothic part, and I could not remember where it was. And I said to myself that would be quite long to do all the street of the neighborhood. And I discovered that it was actually an old problem. In 1962, Meigu Guan was interested in a postman delivering mail to a number of streets such that the total distance walked by the postman was as short as possible. How could the postman ensure that the distance walked was a minimum?

A very close notion is the concept of traversable graph, which is one that can be drawn without taking a pen from the paper and without retracing the same edge. In such a case the graph is said to have an Eulerian trail (yes, from Euler’s bridges problem). An Eulerian trail uses all the edges of a graph. For a graph to be Eulerian all the vertices must be of even order.

An algorithm for finding an optimal Chinese postman route is:

1. List all odd vertices.
2. List all possible pairings of odd vertices.
3. For each pairing find the edges that connect the vertices with the minimum weight.
4. Find the pairings such that the sum of the weights is minimised.
5. On the original graph add the edges that have been found in Step 4.
6. The length of an optimal Chinese postman route is the sum of all the edges added to the total found in Step 4.
7. A route corresponding to this minimum weight can then be easily found.

For the first steps, we can use the codes from Hurley & Oldford’s Eulerian tour algorithms for data visualization and the PairViz package. First, we have to load some R packages

require(igraph) require(graph) require(eulerian) require(GA)

Then use the following function from stackoverflow,

make_eulerian = function(graph){ info = c("broken" = FALSE, "Added" = 0, "Successfull" = TRUE) is.even = function(x){ x %% 2 == 0 } search.for.even.neighbor = !is.even(sum(!is.even(degree(graph)))) for(i in V(graph)){ set.j = NULL uneven.neighbors = !is.even(degree(graph, neighbors(graph,i))) if(!is.even(degree(graph,i))){ if(sum(uneven.neighbors) == 0){ if(sum(!is.even(degree(graph))) &gt; 0){ info["Broken"] = TRUE uneven.candidates &lt;- !is.even(degree(graph, V(graph))) if(sum(uneven.candidates) != 0){ set.j &lt;- V(graph)[uneven.candidates][] }else{ info["Successfull"] &lt;- FALSE } } }else{ set.j &lt;- neighbors(graph, i)[uneven.neighbors][] } }else if(search.for.even.neighbor == TRUE &amp; is.null(set.j)){ info["Added"] &lt;- info["Added"] + 1 set.j &lt;- neighbors(graph, i)[ !uneven.neighbors ][] if(!is.null(set.j)){search.for.even.neighbor &lt;- FALSE} } if(!is.null(set.j)){ if(i != set.j){ graph &lt;- add_edges(graph, edges=c(i, set.j)) info["Added"] &lt;- info["Added"] + 1 } } } (list("graph" = graph, "info" = info))}

Then, consider some network, with 12 nodes

g1 = graph(c(1,2, 1,3, 2,4, 2,5, 1,5, 3,5, 4,7, 5,7, 5,8, 3,6, 6,8, 6,9, 9,11, 8,11, 8,10, 8,12, 7,10, 10,12, 11,12), directed = FALSE)

To plot that network, use

V(g1)$name=LETTERS[1:12] V(g1)$color=rgb(0,0,1,.4) ly=layout.kamada.kawai(g1) plot(g1,vertex.color=V(newg)$color,layout=ly) Then we convert it to some traversable graph by adding 5 vertices eulerian = make_eulerian(g1) eulerian$info broken Added Successfull 0 5 1 g = eulerian$graph as shown below ly=layout.kamada.kawai(g) plot(g,vertex.color=V(newg)$color,layout=ly) We cut those 5 vertices in two part, and therefore, we add 5 artificial nodes

A=as.matrix(as_adj(g)) A1=as.matrix(as_adj(g1)) newA=lower.tri(A, diag = FALSE)*A1+upper.tri(A, diag = FALSE)*A for(i in 1:sum(newA==2)) newA = cbind(newA,0) for(i in 1:sum(newA==2)) newA = rbind(newA,0) s=nrow(A) for(i in 1:nrow(A)){ Aj=which(newA[i,]==2) if(!is.null(Aj)){ for(j in Aj){ newA[i,s+1]=newA[s+1,i]=1 newA[j,s+1]=newA[s+1,j]=1 newA[i,j]=1 s=s+1 }}}

We get the following graph, where all nodes have an even number of vertices !

newg=graph_from_adjacency_matrix(newA) newg=as.undirected(newg) V(newg)$name=LETTERS[1:17] V(newg)$color=c(rep(rgb(0,0,1,.4),12),rep(rgb(1,0,0,.4),5)) ly2=ly transl=cbind(c(0,0,0,.2,0),c(.2,-.2,-.2,0,-.2)) for(i in 13:17){ j=which(newA[i,]&gt;0) lc=ly[j,] ly2=rbind(ly2,apply(lc,2,mean)+transl[i-12,]) } plot(newg,layout=ly2) Our network is now the following (new nodes are small because actually, they don’t really matter, it’s just for computational reasons)

plot(newg,vertex.color=V(newg)$color,layout=ly2, vertex.size=c(rep(20,12),rep(0,5)), vertex.label.cex=c(rep(1,12),rep(.1,5))) Now we can get the optimal path n &lt;- LETTERS[1:nrow(newA)] g_2 &lt;- new("graphNEL",nodes=n) for(i in 1:nrow(newA)){ for(j in which(newA[i,]&gt;0)){ g_2 &lt;- addEdge(n[i],n[j],g_2,1) }} etour(g_2,weighted=FALSE)  "A" "B" "D" "G" "E" "A" "C" "E" "H" "F" "I" "K" "H" "J" "G" "P" "J" "L" "K" "Q" "L" "H" "O" "F" "C"  "N" "E" "B" "M" "A" or edg=attr(E(newg), "vnames") ET=etour(g_2,weighted=FALSE) parcours=trajet=rep(NA,length(ET)-1) for(i in 1:length(parcours)){ u=c(ET[i],ET[i+1]) ou=order(u) parcours[i]=paste(u[ou],u[ou],sep="|") trajet[i]=which(edg==parcours[i]) } parcours  "A|B" "B|D" "D|G" "E|G" "A|E" "A|C" "C|E" "E|H" "F|H" "F|I" "I|K" "H|K" "H|J" "G|J" "G|P" "J|P"  "J|L" "K|L" "K|Q" "L|Q" "H|L" "H|O" "F|O" "C|F" "C|N" "E|N" "B|E" "B|M" "A|M" trajet  1 3 8 9 4 2 6 10 11 12 16 15 14 13 26 27 18 19 28 29 17 25 24 7 22 23 5 21 20 Let us try now on a real network of streets. Like Missoula, Montana. I will not try to get the shapefile of the city, I will just try to replicate the photography above. If you look carefully, you will see some problem : 10 and 93 have an odd number of vertices (3 here), so one strategy is to connect them (which explains the grey line). But actually, to be more realistic, we start in 93, and we end in 10. Here is the optimal (shortest) path which goes through all vertices. Now, we are ready for Halloween, to go through all streets in the neighborhood ! # Halloween and candies (a ballot problem) This year, for Halloween, a post on candies (I promise, next year I will write another post on zombies). But I don’t want to focus on the kids problems (last year, we tried to minimize their walking distance to collect as much candies as possible, with part 1 and part 2), I want to discuss my own problems. Because usually, the kids wear their costumes, and they go in the streets, they knock on the doors, while I stay at home. So I’m the one, with a bag full of candies, waiting for kids to knock on our door, and then I give them some candies (if they wear a costume). Consider the following problem. Assume that we start with red candies, and black ones, with . The thing is that no one like those black candies. What could be the probability that for the kids that will get candies after knocking at my door (with for convenience, but we will also consider the more general case where I have to many candies, , later on), the probability to get a red candy is always larger than the probability to have a black candy ? This is somehow related to the popular ballot problem, proposed (and solved) by Whitworth in 1878, but he wrote it only in the fourth edition of Choice and Chance, in 1886 (this is what the legend told us). In 1887, Joseph Bertrand proposed a similar problem, and Désiré André introduced the reflection principle to solve it. The problem is simple : consider an election between two candidates, A (who receives votes) and B (who receives votes). A wins the election (). If the ballots are cast one at a time, what is the probability that A will lead throughout the voting? For those who don’t remember the conclusion, the probability is here quite simple, Observe that some geometry proofs were given, later on, by Aebly or Mirimanoff, both in 1923, as well as Howard Grossman in the 1950’s (see the discussion on http://academiclogbook.blogspot.ca/…). Actually, http://futilitycloset.com// produced the following geometric proof (with no clear reference), We start at O, where no votes have been cast. Each vote for A moves us one point east and each vote for B moves us one point north until we arrive at E, the final count, (mn). If A is to lead throughout the contest, then our path must steer consistently east of the diagonal line OD, which represents a tie score. Any path that starts by going north, through (0,1), must cut OD on its way to E. If any path does touch OD, let it be at C. The group of such paths can be paired off as p and q, reflections of each other in the line OD that meet at C and continue on a common track to E. This means that the total number of paths that touch OD is twice the number of paths p that start their journey to E by going north. Now, the first segment of any path might be up to m units east or up to units north, so the proportion of paths that start by going north is n/(m + n), and twice this number is 2n/(m + n). The complementary probability — the probability of a path not touching OD — is (m –n)/(m + n). But let’s try to solve our problem. Let and denote the number of black and red candies, respectively after the th kid git his (or her) candy. Yes, one at a time. Here, and . What we want is Using this formulation, we recognize the ballot problem. Almost. Actually, in the original ballot problem (see Bertrand (1887)), we have to compute the probability that one candidate remains strictly ahead the other one throughout the count. With a strict condition, we get the well-known probability (given previously) Here, ties are allowed, and we can prove (easily) that (again, there is some nice geometric interpenetration of that result). It is also possible to get numerically that value using the following function, which will generate a trajectory, and return some indicators (with or without ties) > red_black=function(sd){ + set.seed(sd) + vectcandy=sample(c(rep("R",r),rep("B",b))) + v1=rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B")) + v2=cumsum(rev(vectcandy)=="R")<= cumsum(rev(vectcandy)=="B") + return(list(evol=cbind(rev(cumsum(rev(vectcandy)=="R")), + rev(cumsum(rev(vectcandy)=="B")),v1),list=vectcandy,test=(sum(v1)==0), + ballot=(sum(v2)==0),when=min(which(v1==1))))} (here I compute the ballot-type index, where ties are not allowed, and the candy-type index). If we generate 100,000 scenarios, starting with 50 red and 25 black candies, we get > r=50 > b=25 > M=sapply(1:100000,red_black) ­­ > mean(unlist(M[3,]))  0.50967 which can be compared with the theoretical value > (r+1-b)/(r+1)  0.5098039 We can also get the distribution of the first time we have more black candies than red ones left (given that this event occur) > Z=unlist(M[5,]) > Z=Z[Z<Inf] > hist(Z,breaks=seq(0,80),probability=TRUE,col="light blue", + border=NA,xlab="",main="") There might be some analytically formula that can be derived, but I have to confess that I am becoming extremely lazy, Assume now that this year, kids do not show up at my door (for some reason). Assume that kids show up. We can see how the probability will change, with , > r=50 > b=25 > impact_n = function(n){ + red_black=function(sd,nb=n){ + set.seed(sd) + vectcandy=sample(c(rep("R",r),rep("B",b))) + v=(rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B")))[1:nb] + return(list(list=vectcandy,test=(sum(v)==0),when=min(which(v==1))))} + M=sapply(1:10000,red_black) + return(mean(unlist(M[2,])))} Yes, not only I am too lazy to derive analytic formulas, I am so lazy that I do not try to optimize my code. Here, the evolution of the probability, as a function of is > V=Vectorize(impact_n)(25:75) > plot(25:75,V) Fun isn’t it? But now, I have to conclude my post, to work a little bit on my make-up : I have learnt so many thinks at the Montreal Zombie Walk a few days ago that kids willing to knock at my door will be scared to death. I guess I will keep all the candies for me this year ! # Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 2) Continuons la discussion de ce matin, avec l’exemple des capitales américaines. On peut trouver une telle liste liste en ligne sur wikipedia, > capital=read.table( + "http://freakonometrics.free.fr/US-capital.csv", + header=TRUE,sep=";") > capital=capital[-c(2,11),] On va enlever Hawaï et l’Alaska, qui forcent à faire de trop grands détours. Dans TSP, sous R, il y a des listes de villes américaines, avec les coordonnées. En bricolant un peu, on va récupérer les coordonnées de la plupart des capitales, > library(maps) > library(maptools) > library(sp) > library(TSP) > Lcap=paste(as.character(capital[,1]), + as.character(capital[,2]),sep=", ") > M=as.matrix(USCA312) > L=labels(M)[] > Lcap=Lcap[which(Lcap%in%L)] > k=which(L%in%Lcap) > M=M[k,] > M=M[,k] > listeUSA=TSP(M,L[k]) Et cette fois on est prêt, > tour=solve_TSP(listeUSA, method = "nn") on a lancé l’algorithme, et on peut faire un dessin, > COORD=data.frame(USCA312_coords[k,]) > COORD=COORD[as.numeric(tour),] > map('state',fill=TRUE,col="yellow") > lines(COORD$coords.x1,COORD$coords.x2, + lwd=3,col="blue") > points(COORD$coords.x1,COORD$coords.x2, + pch = 19, cex = 1.2, col = "red") pas mal, n’est pas ? Bon, mais comme d’habitude, rien ne garantie que l’on ait trouvé le trajet optimal… En fait, le soucis est que si on relance l’algorithme, on retombe sur un autre trajet, voire encore un autre si on lancer une nouvelle fois l’algorithme, C’est pénible, n’est ce pas ? En fait, si on lance 1,000 fois l’algorithme on obtient la distribution (des distances optimales) suivante autrement dit, on a encore beaucoup de variabilité, avec une distance optimale qui peut varier d’environ 10% entre deux boucles d’algorithmes. Cela dit, cela reste faible par rapport à un trajet au hasard, Et d’ailleurs, si on compare la distance obtenue en tirant les villes au hasard, on arrive à un taux de l’ordre de 27% (qui peut être comparé à ce que nous avions obtenu en tirant au hasard dans le carré unité) > HAZARD=OPTIMAL=RATIO=rep(NA,500) > for(s in 1:500){ + HAZARD[s]=tsp.longueur(M,sample(1:nrow(M))) + tour=solve_TSP(listeUSA, method = "nn") + OPTIMAL[s]=tsp.longueur(M,as.numeric(tour)) + RATIO[s]=HAZARD[s]/ + OPTIMAL[s] + } Bref, compte tenu de la (grande) dispersion des résultats optimaux obtenu, on pourrait dire que cet algorithme est assez mauvais. Même s’il est très rapide. Un autre avantage est que l’on peut spécifier la ville de départ, par exemple. Ici, si on souhaite partir de Floride, on peut obtenir le trajet suivant > tour=solve_TSP(listeUSA, method = "nn", + control=list(start=41)) Passionnant n’est-ce pas ? Et encore, ça ne va pas m’aider à trouver le parcours optimal pour les enfants, car pour les enfants, il faut spécifier le point de départ (voire le point d’arrivée), il faut en plus rester sur la route (voire traverser aux passages pour piétons). On peut aussi vouloir faire une pause pipi au milieu de la promenade…. Bref, on a pas mal de contraintes ! Mais tout est expliqué dans le livre fabuleux de William Cook, In Pursuit of the Traveling Salesman, qui reprend toute l’histoire des algorithmes les plus fins pour résoudre ce joli problème de maths… # Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 1) Bon, ben cette fois ça y est, c’est l’Halloween… et comme tous les ans, mon rôle devrait se limiter à creuser les citrouilles, décorer la maison… et attendre à la porte de la maison que les enfants des voisins viennent me dévaliser ma réserve de bonbons ! Je pourrais à la rigueur revendiquer un petit quelque chose dans le déguisement de mon fils (on a passé pas mal de temps à feuilleter les premiers tomes de Walking Deads, histoire de mieux saisir l’essence des personnages de morts-vivants, en plus de la Zombie Walk la fin de semaine passée, en centre-ville). Mais cette année, pas de billet sur les zombies comme l’an passé, mais un vrai problème pratique: comment optimiser son trajet, quand on sait où sont les bonnes maisons, où on a de fortes chances d’avoir plein de bonbons (en particulier chez les parents des amis) ? C’est assez proche de l’idée de passer par un maximum de maisons, et ce à une heure raisonnable (histoire qu’il reste des bonbons à récupérer). Pour tous ceux qui ont fait un peu de programmation dynamique, c’est un problème classique, connu sous le nom du “voyageur de commerce“. Ce problème est né en juillet 1954 avec “the shortest route for a traveling salesman” paru dans Newsweek ci-dessous (que l’on retrouve en détails dans le livre 50 Years of Integer Programming) Enfin, disons qu’il a été popularisé à cette époque, car Euler évoquait déjà un problème similaire quand il travaillait sur le problème des sept ponts de Königsberg, dans “Solutio problematis ad geometriam situs pertinentis“, (même si ce problème consiste à étudier l’existence d’un trajet sous des contraintes de passages par des nœuds, et pas vraiment de l’optimisation). On retrouve aussi un exemple en 1832, en Allemagne, dans un livre intitulé Der Handlungsreisende – wie er sein soll und was er zu thun hat, um Auftraege zu erhalten und eines gluecklichen Erfolgs in seinen Geschaeften gewiss zu sein – Von einem alten Commis-Voyageur où un trajet passant par 45 villes en Allemagne et en Suisse devait être optimisé (le schéma ci-dessus montre celle décrite en 1832, pour un total de 1,285km, mais il est possible de faire une boucle de 1,248km). En fait le problème a explosé en 1962 (à cette époque, il semble que les gens pouvaient encore penser que faire des maths pouvait être le fun), lors que Procter & Gamble ont offert un prix de 10,000$ à celui qui trouverait le chemin le plus court pour que Toody et Muldoom, les conducteurs de la voiture 54 dans une série populaire à l’époque, passent par 33 villes, aux États-Unis. Cela dit, Carl Menger avait noté dès 1930 a noté que ce problème devrait etre difficile à résoudre. Facile à écrire, mais difficile à résoudre explicitement. Dans son livre, William Cooke rappelle qu’avec 33 villes, il y aurait

131,565,418,466,846,765,083,609,006,080,000,000

chemins possibles, dont il faudrait calculer la longueur. Et qu’avec l’ordinateur le plus puissant en 2009 (effectuant 1.5 million de milliards d’opérations à la seconde) il faudrait 28,000 milliards d’année pour mener à bien les calculs. Même si on essayait de tenir compte des progrès de l’informatique dans les mille prochaines années, ça risque de prendre du temps !..

En 1967, Jack Edmonds notait “I conjecture that there is no good algorithm for the traveling salesman problem” où “good” est une manière élégante pour dire un temps polynomial en , le nombre de villes. Ce qui s’appelle les problèmes de la classe . Il semblerait que le problème du voyageur de commerce soit un problème  et pas de la classe (cf le chapitre du livre de Sanjeev Arora et Boaz Barak, ou le chapitre 11 des notes de cours de Olivier Bournez, que j’ai pu découvrir grâce à @dmonniaux, ou formellement de la classe  et pas , i.e.  complet), voir aussi la discussion sur wikipedia),

Et depuis, l’institut Clay a proposé $1,000,000 à celui qui prouverait que les problèmes ne sont pas de type (ou qui prouveraient qu’ils le sont… laconjecture reste ouverte, semble-t-il). Bref, les enjeux financiers derrière ont considérablement augmenté (on pourra consulter Sipser (2007) pour la petite histoire) Bon, et si l’histoire est intéressante, ce n’est pas une raison pour la lire en se croisant les bras… essayons de programmer un algorithme (simple) afin d’approcher une solution, car je ne pourrais me contenter de dire à mes enfants “c’est un problème complet, papa a capitulé“. Commençons par tirer 5 points au hasard dans le carré unité, > n = 5 > (x = matrix(runif(2*n),nr=n)) [,1] [,2] [1,] 0.4142621 0.35600725 [2,] 0.9477231 0.66938653 [3,] 0.2518656 0.19900873 [4,] 0.5314231 0.22954670 [5,] 0.5156522 0.09066741 et commençons par calculer la distance entre les points (i.e. une matrice, ou ici une matrice triangulaire inférieure) > (d = dist(x)) 1 2 3 4 2 0.6186980 3 0.2258786 0.8399243 4 0.1723919 0.6056111 0.2812205 5 0.2840514 0.7222195 0.2851688 0.1397719 On va mettre ce triangle sous forme de matrice (pleine), ça sera plus simple à manipuler par la suite > (d = as.matrix(d)) 1 2 3 4 5 1 0.0000000 0.6186980 0.2258786 0.1723919 0.2840514 2 0.6186980 0.0000000 0.8399243 0.6056111 0.7222195 3 0.2258786 0.8399243 0.0000000 0.2812205 0.2851688 4 0.1723919 0.6056111 0.2812205 0.0000000 0.1397719 5 0.2840514 0.7222195 0.2851688 0.1397719 0.0000000 Le principe, pour initialiser l’algorithme est de tirer au hasard un trajet, c’est à dire une permutation de nos points, > (o = sample(1:n))  2 4 1 5 3 Il faut alors calculer le 4 distances entre ces 5 villes (on ne prévoit pas de retour à la ville de départ, pour l’instant). Le plus simple est de jouer avec cette matrice de la réarranger, > d[o[1:(n-1)],][,o[2:n]] 4 1 5 3 2 0.6056111 0.6186980 0.7222195 0.8399243 4 0.0000000 0.1723919 0.1397719 0.2812205 1 0.1723919 0.0000000 0.2840514 0.2258786 5 0.1397719 0.2840514 0.0000000 0.2851688 et de prendre les 4 valeurs sur la diagonale > diag(d[o[1:(n-1)],][,o[2:n]])  0.6056111 0.1723919 0.2840514 0.2851688 Ce sont les 4 distances que l’on cherche, entre les 5 villes. La distance totale parcourue est alors > sum(diag(d[o[1:(n-1)],][,o[2:n]]))  1.347223 L’idée est ensuite de permuter deux des villes > (i=sample(1:n,2))  3 1 (ici la première et la troisième de la liste, mais on pourrait décider de fixer la première si on doit partir de là) et de voir si la nouvelle distance est plus courte, ou pas, > os=o > os[i]=o[rev(i)] > sum(diag(d[os[1:(n-1)],][,os[2:n]]))  1.785391 La distance étant plus longue, on ne va pas permuter. On pourrait tenter 1000 fois de permuter deux villes. Si on ne parvient pas à réduire la distance totale, c’est qu’on a atteint une valeur intéressante (à défaut de pouvoir prouver qu’elle est optimale). Et si on peut améliorer la distance, on permute. Le plus simple est donc de faire des petites fonctions. La première reprend le calcul de la distance totale, sur la diagonale de la matrice permutée, à matrice de distance donnée et à permutation donnée, > tsp.longueur=function(matrice,ordres) { + n=length(ordres) + sum(diag(matrice[ordres[1:(n-1)],][,ordres[2:n]])) } La seconde est juste une boucle, > tsp.optimal=function(d,N=1000) { + d=as.matrix(d) + n=ncol(d) + o=sample(1:n) + v=tsp.longueur(d,o) + k=0 + while(k < N) { + i<-sample(1:n,2) + os=o + os[i]=o[rev(i)] + w=tsp.longueur(d,os) + if(w < v) { + v=w + o=os + k=0} else {k=k+1}} + list(ordre=o,longueur=v)} On peut se lancer maintenant, avec par exemple 30 points simulés, sur le carré unité, > set.seed(1) > n=30 > x=matrix(runif(2*n),nr=n) > o=sample(1:n) > plot(x,xlim=0:1,ylim=0:1,xlab="",ylab="") > lines(x[o,],col='blue') > r=tsp.optimal(dist(x)) > os=r$ordre
>  lines(x[os,],col='blue',lwd=1.5)
avec la première trajectoire, à gauche, qui reste en trait fin sur le dessin de droite,  Pas mal, non ? Bon, en fait le soucis est que si on joue plusieurs fois avec notre fonction, on obtient toujours un graphique différent,  Autrement dit, je trouve (certes) des chemins pas trop long, mais je suis loin de trouver le chemin optimal. Et histoire de clotûrer le débat, notons qu’on peut faire un peu de simulation, afin de comparer la distance obtenue en faisant le trajet au hasard, et un trajet “optimal“,

> RATIO=HASARD=OPTIMAL=matrix(NA,500,10)
> for(m in 1:10){
+ n=5*m
+ for(s in 1:500){
+  x=matrix(runif(2*n),nr=n)
+  r=tsp.opt(dist(x))
+  HAZARD[s,m]=tsp.longueur(as.matrix(dist(x)),1:n)
+  OPTIMAL[s,m]=r\$longueur
+  RATIO[s,m]=HAZARD[s,m]/
+             OPTIMAL[s,m]
+ }
+}

Si on compare la longueur moyenne du trajet (sur 500 scénarios) avec un tour au hasard (en rouge) et ou tour optimisé (en bleu), on a ce qui donne un ratio moyen (ou une moyenne des ratios, les deux sont représentés) Autrement dit, avec une dizaine de points éparpillés (au hasard) dans le carré unité, un tour au hasard ne sera que deux fois plus long… étonnant non ?

Mais on doit pouvoir aller plus loin, parce que sous R, il y a des fonctions (et des packages) pour faire des algorithmes plus compliqués… à suivre avant la tournée des bonbons de demain soir !

# Halloween: zombies walk versus random walk

I don’t know if it’s because of Halloween (or the new season of walking dead), but I had weird thoughts recently. For instance, I wondered if zombies walk was a random walk… Note that I am not the only one: there are papers on that topic on Wired or even Nature (at least on one of their blogs). Before modeling that walk, it might be better to describe the process, and to generate zombies walk through simulations.