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 inPython 3.7.3, which mostly means it hastyping, 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.

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

- Find the center point of all the places I want to visit.
- Sort the places I want to visit according to their angle from the center.
- 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:

classGeoPoint:

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:

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

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

defget_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

defvisualize(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:

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

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

defplot_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