Creatie door Martin Fiers (mfiers@gmail.com). Vrij om te gebruiken voor alle doeleinden. Bewust in het Nederlands geschreven met Nerdland podcast luisteraars in het achterhoofd (ik kan Engels hoor ;) ).
Sam Gooris gaat in zijn liedje 'Ambiance Ambiance' heel veel steden aandoen (24 om precies te zijn). Echter, de volgorde waarin deze worden aangedaan is ver van optimaal. Dit werd opgemerkt door Jeroen Baert in de nerdland podcast (maandoverzicht maart 2019), en er werd een warme oproep gedaan om een goeie volgorde te vinden.
Bij deze...
https://www.songteksten.nl/songteksten/86757/sam-gooris/ambiance%2C-ambiance%2C-ambiance.htm
Naar een bal in Mal
En naar een tent in Gent
We vieren feest in Leest
En maken sfeer in Peer
We geven gas in As
En we slaan tilt in Tielt
We worden zot in Lot
En tureluurs in Puurs, komaan
Ambiance, ambiance, ambiance
Als het weekend is gaan we uit de bol
Ambiance, ambiance, ambiance
Met de kliek op zwier en we maken lol
Een frisse pint in Lint
Kip met rijst in Heist
Dan naar een keet in Reet
En op cafe in Bree
We drinken kriek in Schriek
We kijken scheel in Geel
Dan wordt het leuk in Leut
En wordt het zwoel in Doel
Ambiance, ambiance, ambiance
Als het weekend is gaan we uit de bol
Ambiance, ambiance, ambiance
Met de kliek op zwier en we maken lol, hey
Een dikke knuffel in Duffel
Een aai in Sinnaai
Een mooie borst in Vorst
En sex-appeal in Niel
Uit de klere in Bere
Uit de rits in Gits
Een natte droom in Boom
Een blonde maagd in Haacht
Ambiance, ambiance, ambiance
Als het weekend is gaan we uit de bol
Ambiance, ambiance, ambiance
Met de kliek op zwier en we maken lol
Tof liedje ;). Eerst encoderen we de steden, de locaties (latitude, longitude) zijn eenvoudig gevonden op https://www.latlong.net/. Voor 24 steden wou ik nog net niet de moeite doen om code te schrijven om via een API al het opzoekwerk te doen en google accounts aan te maken, dus heb ik maar manueel alles opgezocht:
cities = [
('Mal', (50.768450, 5.521440)),
('Gent', (51.054340, 3.717424)),
('Leest', (51.033699, 4.419630)),
('Peer', (51.185883, 5.289227)),
('As', (51.000542, 5.572203)),
('Tielt', (51.000061, 3.324990)),
('Lot', (50.765170, 4.276202)),
('Puurs', (51.071998, 4.286193)),
('Lint', (51.127670, 4.495450)),
('Heist', (51.337436, 3.239979)), # Aan zee.
('Reet', (51.104485, 4.405890)),
('Bree', (51.141090, 5.597620)),
('Schriek', (51.021873, 4.709033)),
('Geel', (51.162319, 4.990880)),
('Leut', (50.989192, 5.741560)),
('Doel', (51.317621, 4.245205)),
('Duffel', (51.095718, 4.506199)),
('Sinnaai', (50.780688, 4.792712)),
('Vorst', (50.809143, 4.317751)),
('Niel', (51.109865, 4.330340)),
('Bere', (50.949909, 3.288330)), # Meulebeke
('Gits', (50.997592, 3.110329)),
('Boom', (51.087379, 4.366722)),
('Haacht', (50.976929, 4.638682)),
]
import numpy as np
city_list = np.array([(x, y) for _, (x, y) in cities])
distance = 0
nr_cities = len(cities)
route = [x for x in range(nr_cities)]
for i in range(nr_cities):
distance += np.sqrt((city_list[route[i]][0] - city_list[route[(i + 1) % nr_cities]][0] ) ** 2 + (city_list[route[i]][1] - city_list[route[(i + 1) % nr_cities]][1] ) ** 2)
km_per_angle = 111
print(f"Aantal steden: {len(cities)}. Afstand in het originele liedje: {distance * km_per_angle} km.")
Aantal steden: 24. Afstand in het originele liedje: 2210.2817842249287 km.
2210 km. Dat is veel. Eventjes zien hoe Sam zijn traject heeft afgelegd:
from ipyleaflet import Marker, Map, Polyline
center = list(np.mean(city_list, axis=0))
m = Map(center=center, zoom=9)
for city, center in cities:
marker = Marker(location=center, title=city, draggable=False)
m.add_layer(marker)
line = Polyline(locations=[[
[[city_list[x][0], city_list[x][1]]
for x in range(len(cities))]
]])
m.add_layer(line)
m
Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'max_zoom': 19, 'attribution': 'Map …
Dat kan beter!
Ik heb aan 2 algoritmes gedacht: simulated annealing en een genetische algoritme. Ik heb voor simulated annealing gekozen omdat het toch net iets sneller te implementeren is dan een genetisch algoritme (dat laatste had nog cooler geweest, dan had ik ook Hetty's aandacht getrokken ;)).
(voor implementatiedetails zie vanonder deze notebook)
%pylab inline
import numpy as np
import random
import numba as nb
Populating the interactive namespace from numpy and matplotlib
@nb.njit()
def mix_route(route, nr_cities):
"""Generate a variation of the route (in-place).
Very naive implementation that just swaps two items in the route"""
i, j = np.random.choice(np.arange(0, nr_cities), 2, replace=False)
tmp = route[j]
route[j] = route[i]
route[i] = tmp
@nb.njit()
def simulated_annealing_tsp(samples, city_list, start_t=10000., end_t=1.):
nr_cities = len(city_list)
# Choose random route to start with
route = np.random.choice(np.arange(0, nr_cities), nr_cities, replace=False)
start_route = route
distance_history = np.zeros((samples,))
prob_change = np.zeros((samples,))
best_distance = 1e100
decay_rate = 10 ** ((np.log10(end_t) - np.log10(start_t)) / (samples - 1))
t = start_t
for idx in range(samples):
new_route = np.copy(route)
mix_route(new_route, nr_cities)
# Numba-enabled
old_dist = 0
new_dist = 0
for i in range(nr_cities):
old_dist += np.sqrt((city_list[route[i]][0] - city_list[route[(i + 1) % nr_cities]][0] ) ** 2 + (city_list[route[i]][1] - city_list[route[(i + 1) % nr_cities]][1] ) ** 2)
new_dist += np.sqrt((city_list[new_route[i]][0] - city_list[new_route[(i + 1) % nr_cities]][0] ) ** 2 + (city_list[new_route[i]][1] - city_list[new_route[(i + 1) % nr_cities]][1] ) ** 2)
prob = np.exp((old_dist - new_dist) / t)
if prob > random.random():
# accept new distance
route = new_route
distance_history[idx] = old_dist
prob_change[idx] = prob
t = t * decay_rate
return start_route, distance_history, prob_change, route
import time
t_start = time.time()
start_route, distance_history, prob_change, route = simulated_annealing_tsp(samples=10000000,
start_t=100.,
end_t=0.01,
city_list=city_list)
t_end = time.time()
print(f"Simulation time: {t_end - t_start}")
Simulation time: 16.145136833190918
print(f"Beste afstand na {len(distance_history)} iteraties: {distance_history[-1] * km_per_angle} km.")
Beste afstand na 10000000 iteraties: 703.0181067123403 km.
Om te zien wat het algoritme doet (beetje), kunnen we eens zien hoe de afstand evolueert tijdens het annealen:
title("Afstand i.f.v. iteraties")
xlabel("Iteratie")
ylabel("Afstand (km)")
_ = plot(distance_history[::100] * km_per_angle)
Eens kijken hoe het traject er nu uit ziet:
from ipyleaflet import Marker, Map, Polyline
center = list(np.mean(city_list, axis=0))
m = Map(center=center, zoom=9)
for city, center in cities:
marker = Marker(location=center, title=city, draggable=False)
m.add_layer(marker)
line = Polyline(locations=[[
[[city_list[route[x]][0], city_list[route[x]][1]]
for x in range(len(cities))]
]])
m.add_layer(line)
m
Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'max_zoom': 19, 'attribution': 'Map …
Met dit algoritme is het doenbaar om een route te vinden rond de 700 km (in vogelvlucht wel te verstaan, maar dat waren ook de regels ;) ). 10 miljoen iteraties is doenbaar op een laptop en daar kan je ook al 50-100 steden mee oplossen.
Als bijkomstig interessant weetje valt het sterk op dat de meeste steden zich rond Antwerpen hebben geconcentreerd. Sam is blijkbaar een Antwerpenaar puur op basis van deze data.
Het simulated annealing algoritme is erg eenvoudig om te implementeren. Inspiratie gehaald van https://ericphanson.com/blog/2016/the-traveling-salesman-and-10-lines-of-python/.
Daarnaast heb ik numba gebruikt omdat het de code te accellereren. Numba is een just-in-time compiler die Python code omzet naar llvm (low-level virtual machine), en dan gecompileerd naar machine code. Geweldige tool. Voor 100000 iteraties duurt het in pure Python 47s, met numba 0.136 seconden (factor 345!). Ik kon dus makkelijk 10 miljoen iteraties doen in ~ 12 seconden.
Daarnaast is het creeeren van nieuwe routes tijdens het algoritme zeer naief. i.p.v. twee routes te wisselen zou je iets meer variatie kunnen toelaten, om sneller te 'scannen' naar verschillende plaatsen in de oplossingsruimte.
Er zijn natuurlijk veel professionelere tools / algoritmes, dus gebruik dit zelf niet als je zelf een TSP wilt oplossen (uiteraard).
Voila.
Martin
mfiers@gmail.com