Today I want to tell you about yet another solution I came up with for the famous Euclidean Traveling Salesperson Problem (E-TSP). The advantages of this heuristic solution is that it is intuitive and thus simple to implement.

There are better (more exact + faster) methods out there, but I chose to take the time and share with you my simple solution since I think it can serve as a great intro to this field, and also as a fresh point of view for my peers in the Smart Mobility domain, who I believe, like me, sometimes just want to get things done.

If you just want to see the code, because, hey, you too just want to get things done, you are welcome to go directly here.

Note: The code in this post is written in Python 3.7.3, which mostly means it has typing, which is my new favourite Python feature 💜

Let’s Get Things Done

Here is a short list of stuff I wanted to do on a Friday morning before my sleepy hometown goes dead for the Israeli weekend: visit the Persian couple who make amazing Tahdig, eat delicious pita from the lovely Druze family, buy some food and medicine for myself and my dog, and also rest for a while in the dog park.

And did I mention already this is a sleepy town? So I have to find the fastest route through all my errands, before all the shops close.

Since E-TSP is an NP-hard problem, I want to implement a heuristic method instead of looking for the optimal solution (yeah, I know this input is of size 7, we’re gonna explore larger inputs at the end).

Done Is Better Than Perfect

It’s true that done is better than perfect, but how can we get this one done??

In the last couple of years I am working as a Smart Mobility Algorithm Developer, and so I tend to get lost sometimes in an endless crawl over mobility algorithms paper, and this time I bumped into this one which suggested that humans are pretty good and finding solution to E-TSP with their meer visual intuition, which got me thinking that I can formulize this intuition into code.

source: https://xkcd.com/2173/

So here is the outline of the solution we are going to implement:

  1. Find the center point of all the places I want to visit.
  2. Sort the places I want to visit according to their angle from the center.
  3. Perform a greedy local search to optimize this angular route.

Ok so Let’s Get Things Done Already

1. Laying Foundations

We will use the following packages:

import numpy as np
import math

and define the following GeoPoint class we will be working with:

class GeoPoint:
def __init__(self, lng: float, lat: float):
# Why 5 digits? It's 1 m. accuracy according to
# https://en.wikipedia.org/wiki/Decimal_degrees
self.lng = round(lng, 5)
self.lat = round(lat, 5)

def __repr__(self):
# Copy-pastable format for most map applications
return f"[{self.lng}, {self.lat}]"

And from now on we will assume that the list of places we want to visit are held in geo_points: [GeoPoint] (in Python 3.7.3 this means that geo_points is a list of instances of GeoPoint).

2. Finding the Angular Route

We first want to find the center:

def get_geo_point_of_center(geo_points: [GeoPoint]) -> GeoPoint:
lng_list, lat_list = list(zip(*[[g.lng, g.lat]
for g in geo_points]))
lng_center, lat_center = np.mean(np.array([lng_list, lat_list]),
axis=1)
return GeoPoint(lng_center, lat_center)

and then compute the angle from the center:

def get_angle_from_reference_geo_point_in_deg(
reference_geo_point: GeoPoint,
other_geo_point: GeoPoint) -> float:
    x = other_geo_point.lng - reference_geo_point.lng
y = other_geo_point.lat - reference_geo_point.lat
angle_from_reference_in_deg = math.degrees(math.atan2(y, x))
    return angle_from_reference_in_deg

Finally, wrap it up using the built-in sorted function to get the angular route:

def get_angular_route_idxs(geo_points: [GeoPoint]) -> [int]:
center = get_geo_point_of_center(geo_points)
route_idxs = sorted(
list(range(len(geo_points))),
key=lambda i:
get_angle_from_reference_geo_point_in_deg(
center, geo_points[i]),
reverse=True)
return route_idxs

Note: we want to keep the indices of the angular route.

2.5. Some Visualization Tips

If you don’t want to manually draw your route on top of an s2map, like I “lazily” did above, here are some basic tips for using the wonderful geojson format, and the even more wonderful geojson.io visualization service and accompanying package.

In order to not make this into a geojson tutorial, I will just spill it out in its final form, so you can copy-paste this happily:

from geojson import LineString, Point, Feature, FeatureCollection
def visualize(route_idxs: [int], geo_points: [GeoPoint]):
lng_lat_list = [tuple([geo_points[i].lng, geo_points[i].lat])
for i in route_idxs]
route = Feature(geometry=LineString(lng_lat_list),
properties={"name": "This is our route",
"stroke": "#8B0000"})
places = [Feature(geometry=Point(lng_lat),
properties={"name": f"Place {route_idxs[i]}",
"symbol": str(i)[-1],
"marker-color": "#00008B"})
for i, lng_lat in enumerate(lng_lat_list)]
    features = [route] + places
feature_collection = FeatureCollection(features)

return feature_collection

And if you are working with jupyter notebook you can use these lines to automagically open a geojson.io tab with your visualization:

import geojsonio
import json
route_idxs = get_angular_route_idxs(geo_points)
feature_collection = visualize(route_idxs, geo_points)
geojsonio.display(json.dumps(feature_collection));

Expected result:

3. Optimizing

In our greedy local search we will:
- Perform a fixed number of iterations.
- In every iteration we will:
* Change the index of one random place in the list,
* Keep this change if the new route is shorter.

To be able to do this, we will assume we hold a distances_array which is a numpy array holding in location i,j the distance between places i and j .

Using this array we can compute the length of a given route:

def get_route_len(distances_array: np.array, route_idxs):
route_len = sum([distances_array[i1][i2]
for i1, i2 in zip(route_idxs[:-1],
route_idxs[1:])])
return route_len

Finally our simple optimization looks like this:

import random
random_seed = 42
random.seed(random_seed)
def optimize_route(distances_array: np.array, 
route_idxs: [int], n_iter: int) -> [int]:
prev_cost = get_route_len(distances_array, route_idxs)

all_idxs = list(range(len(route_idxs)))
for _ in range(n_iter):
i1, i2 = random.sample(all_idxs, 2)
route_idxs[i2], route_idxs[i1] = (route_idxs[i1],
route_idxs[i2])
new_cost = get_route_len(distances_array, route_idxs)
if new_cost < prev_cost:
prev_cost = new_cost
else:
route_idxs[i2], route_idxs[i1] = (route_idxs[i1],
route_idxs[i2])
return route_idxs

And to wrap everything up we can define:

def plot_best_route(geo_points: [GeoPoint], 
distances_array: np.array, n_iter: int) -> None:
route_idxs = get_angular_route(geo_points)
route_idxs = optimize_route(distances_array, route_idxs, n_iter)
feature_collection = visualize(route_idxs, geo_points)
geojsonio.display(json.dumps(feature_collection));

Expected result:

Would you like to see all the code in one place? You are welcome to go here.

Are We Done?

In the jupyter notebook I prepared, with the full pipeline from this post, you can also find a toy example for a larger input, and see that also for 50 points, which is not feasible to solve to optimality, you can get in less than 1 second, a not-so-bad-looking outcome like this one:

Hope you enjoyed taking this journey with me, and stay tuned for more algorithmic posts coming up in the next couple of weeks :D