Genetic algorithm for Traveling Salesman

I have implemented a genetic algorithm in python 3 for a programming assignment, and I think all the logic is correct. A friend of mine has also implemented one which carries out similar logic, however his was done in Java. It took his around 5 seconds to complete 5000 iterations, whereas mine is taking nearly four minutes!

I have tried to determine the bottleneck using the Python profiler (as described in this answer), but all of them seem to be taking quite long:

Are there any optimisations I can make which will improve the duration of my program? Here is the code.

import sys, math, random, heapq
import matplotlib.pyplot as plt
from itertools import chain

if sys.version_info < (3, 0):
sys.exit("""Sorry, requires Python 3.x, not Python 2.x.""")

class Graph:

def __init__(self, vertices):
self.vertices = vertices
self.n = len(vertices)

def x(self, v):
return self.vertices[v][0]

def y(self, v):
return self.vertices[v][1]

# Lookup table for distances
_d_lookup = {}

def d(self, u, v):
"""Euclidean Metric d_2((x1, y1), (x2, y2))"""

# Check if the distance was computed before
if (u, v) in self._d_lookup:
return self._d_lookup[(u, v)]

# Otherwise compute it
_distance = math.sqrt((u[0] - v[0])**2 + (u[1] - v[1])**2)

self._d_lookup[(u, v)], self._d_lookup[(v, u)] = _distance, _distance
return _distance

def plot(self, tour=None):
"""Plots the cities and superimposes given tour"""

if tour is None:
tour = Tour(self, [])

_vertices = [self.vertices[0]]

for i in tour.vertices:
_vertices.append(self.vertices[i])

_vertices.append(self.vertices[0])

plt.title("Cost = " + str(tour.cost()))
plt.plot(*zip(*_vertices), '-r')
plt.scatter(*zip(*self.vertices), c="b", s=10, marker="s")
plt.show()

class Tour:

def __init__(self, g, vertices = None):
"""Generate random tour in given graph g"""

self.g = g

if vertices is None:
self.vertices = list(range(1, g.n))
random.shuffle(self.vertices)
else:
self.vertices = vertices

self.__cost = None

def cost(self):

if self.__cost is None:
self.__cost = 0
for i, j in zip([0] + self.vertices, self.vertices + [0]):
self.__cost += self.g.d(self.g.vertices[i], self.g.vertices[j])
return self.__cost

class GeneticAlgorithm:

def __init__(self, g, population_size, k=5, elite_mating_rate=0.5,
mutation_rate=0.015, mutation_swap_rate=0.2):
"""Initialises algorithm parameters"""

self.g = g

self.population = []
for _ in range(population_size):
self.population.append(Tour(g))

self.population_size = population_size
self.k = k
self.elite_mating_rate = elite_mating_rate
self.mutation_rate = mutation_rate
self.mutation_swap_rate = mutation_swap_rate

"""Implements ordered crossover"""

size = len(mum.vertices)

# Choose random start/end position for crossover
alice, bob = [-1] * size, [-1] * size
start, end = sorted([random.randrange(size) for _ in range(2)])

# Replicate mum's sequence for alice, dad's sequence for bob
for i in range(start, end + 1):
alice[i] = mum.vertices[i]

# Fill the remaining position with the other parents' entries

for i in chain(range(start), range(end + 1, size)):

while mum.vertices[current_mum_position] in bob:
current_mum_position += 1

bob[i] = mum.vertices[current_mum_position]

# Return twins

def mutate(self, tour):
"""Randomly swaps pairs of cities in a given tour according to mutation rate"""

# Decide whether to mutate
if random.random() < self.mutation_rate:

# For each vertex
for i in range(len(tour.vertices)):

# Randomly decide whether to swap
if random.random() < self.mutation_swap_rate:

# Randomly choose other city position
j = random.randrange(len(tour.vertices))

# Swap
tour.vertices[i], tour.vertices[j] = tour.vertices[j], tour.vertices[i]

def select_parent(self, k):
"""Implements k-tournament selection to choose parents"""
tournament = random.sample(self.population, k)
return max(tournament, key=lambda t: t.cost())

def evolve(self):
"""Executes one iteration of the genetic algorithm to obtain a new generation"""

new_population = []

for _ in range(self.population_size):

# K-tournament for parents

# Mate in an elite fashion according to the elitism_rate
if random.random() < self.elite_mating_rate:
if alice.cost() < mum.cost() or alice.cost() < dad.cost():
new_population.append(alice)
if bob.cost() < mum.cost() or bob.cost() < dad.cost():
new_population.append(bob)

else:
self.mutate(alice)
self.mutate(bob)
new_population += [alice, bob]

# Add new population to old
self.population += new_population

# Retain fittest
self.population = heapq.nsmallest(self.population_size, self.population, key=lambda t: t.cost())

def run(self, iterations=5000):
for _ in range(iterations):
self.evolve()

def best(self):
return max(self.population, key=lambda t: t.cost())

# Test on berlin52: http://elib.zib.de/pub/mp-testdata/tsp/tsplib/tsp/berlin52.tsp
g = Graph([(565.0, 575.0), (25.0, 185.0), (345.0, 750.0), (945.0, 685.0),
(845.0, 655.0), (880.0, 660.0), (25.0, 230.0), (525.0, 1000.0),
(580.0, 1175.0), (650.0, 1130.0), (1605.0, 620.0), (1220.0, 580.0),
(1465.0, 200.0), (1530.0, 5.0), (845.0, 680.0), (725.0, 370.0),
(145.0, 665.0), (415.0, 635.0), (510.0, 875.0), (560.0, 365.0),
(300.0, 465.0), (520.0, 585.0), (480.0, 415.0), (835.0, 625.0),
(975.0, 580.0), (1215.0, 245.0), (1320.0, 315.0), (1250.0, 400.0),
(660.0, 180.0), (410.0, 250.0), (420.0, 555.0), (575.0, 665.0),
(1150.0, 1160.0), (700.0, 580.0), (685.0, 595.0), (685.0, 610.0),
(770.0, 610.0), (795.0, 645.0), (720.0, 635.0), (760.0, 650.0),
(475.0, 960.0), (95.0, 260.0), (875.0, 920.0), (700.0, 500.0),
(555.0, 815.0), (830.0, 485.0), (1170.0, 65.0), (830.0, 610.0),
(605.0, 625.0), (595.0, 360.0), (1340.0, 725.0), (1740.0, 245.0)])

ga = GeneticAlgorithm(g, 100)
ga.run()

best_tour = ga.best()
g.plot(best_tour)

• This is an excellent first question, welcome on CodeReview! :) – IEatBagels May 22 '18 at 14:31
• @TopinFrassi Thank you! It gives friendlier vibes than StackOverflow. – Luke Collins May 22 '18 at 14:36

    def x(self, v):
return self.vertices[v][0]

def y(self, v):
return self.vertices[v][1]


I think these two methods are dead code and could be deleted.

    # Lookup table for distances
_d_lookup = {}

def d(self, u, v):
"""Euclidean Metric d_2((x1, y1), (x2, y2))"""

# Check if the distance was computed before
if (u, v) in self._d_lookup:
return self._d_lookup[(u, v)]

# Otherwise compute it
_distance = math.sqrt((u[0] - v[0])**2 + (u[1] - v[1])**2)

self._d_lookup[(u, v)], self._d_lookup[(v, u)] = _distance, _distance
return _distance


Would be good to have better documentation: what's the relationship between u and (x1, y1)?

In my testing, the lookup didn't actually provide any speedup. I think this is mainly because it's such a big key. If I change u and v to be indices into self.vertices, there is a speed saving. Obviously this also means changing Tour.cost, which is the only method which calls Graph.d.

        if vertices is None:
self.vertices = list(range(1, g.n))
random.shuffle(self.vertices)
else:
self.vertices = vertices


I had to reverse engineer from the if case what the meaning of vertices is. A comment explaining it would have been helpful.

    def cost(self):

if self.__cost is None:
self.__cost = 0
for i, j in zip([0] + self.vertices, self.vertices + [0]):
self.__cost += self.g.d(self.g.vertices[i], self.g.vertices[j])


Yakym's proposed code change was buggy (although it's now fixed), but the point about avoiding creating new lists is a good one. An alternative way of fixing Yakym's code which takes into account the goal of avoiding list creation and is further adapted to correspond to my point about using indices as arguments to Graph.d is

            self.__cost = self.g.d(0, self.vertices[0]) + \
sum(map(self.g.d, self.vertices, self.vertices[1:])) + \
self.g.d(self.vertices[-1], 0)


Separate point on cost: since you're always going to call cost() (if nothing else then in the heapq selection), is there any benefit to making it lazy?

    def crossover(self, mum, dad):
"""Implements ordered crossover"""

size = len(mum.vertices)


Why not this?

        size = g.n - 1


        # Choose random start/end position for crossover
alice, bob = [-1] * size, [-1] * size
start, end = sorted([random.randrange(size) for _ in range(2)])

# Replicate mum's sequence for alice, dad's sequence for bob
for i in range(start, end + 1):
alice[i] = mum.vertices[i]

# Fill the remaining position with the other parents' entries

for i in chain(range(start), range(end + 1, size)):

while mum.vertices[current_mum_position] in bob:
current_mum_position += 1

bob[i] = mum.vertices[current_mum_position]


Here we have a familiar culprit when code runs too slow: in list. Testing whether a list contains a value takes linear time: if you want a fast in test, use a set:

        skipalice, skipbob = set(alice), set(bob)
for i in chain(range(start), range(end + 1, size)):

while mum.vertices[current_mum_position] in skipbob:
current_mum_position += 1

bob[i] = mum.vertices[current_mum_position]
current_mum_position += 1


But actually a faster approach seems to be to use a comprehension and advanced indexing:

    def crossover(self, mum, dad):
"""Implements ordered crossover"""

size = g.n - 1

# Choose random start/end position for crossover
start, end = sorted([random.randrange(size) for _ in range(2)])

# Identify the elements from mum's sequence which end up in alice,
# and from dad's which end up in bob
mumxo = set(mum.vertices[start:end+1])

# Take the other elements in their original order
alice = [i for i in dad.vertices if not i in mumxo]
bob = [i for i in mum.vertices if not i in dadxo]

# Insert selected elements of mum's sequence for alice, dad's for bob
alice[start:start] = mum.vertices[start:end+1]

# Return twins


    def mutate(self, tour):


I get a moderate speedup by pulling out n = len(tour.vertices), saving up to $n^2$ calls of len.

    def select_parent(self, k):
"""Implements k-tournament selection to choose parents"""
tournament = random.sample(self.population, k)
return max(tournament, key=lambda t: t.cost())


I haven't tried this, but I wonder whether it might be faster to sort self.population and then sample from range(self.population_size).

        # Retain fittest
self.population = heapq.nsmallest(self.population_size, self.population, key=lambda t: t.cost())


The documentation for heapq.nsmallest says

[It performs] best for smaller values of n. For larger values, it is more efficient to use the sorted() function.

Here n is len(self.population)/2, and I see a small speedup from

        self.population = sorted(self.population, key=lambda t: t.cost())[0:self.population_size]


It might be worth going for a linear quickselect approach instead.

• NB I've been doing performance testing on the smaller test case elib.zib.de/pub/mp-testdata/tsp/tsplib/tsp/burma14.tsp because I was using an online code execution environment with a 60-second time limit; but on that test case the changes mentioned above got the runtime down from 40 seconds to 28 seconds. – Peter Taylor May 23 '18 at 15:56
• Thank you very much for all these! I'll be sure to read through them soon and give you feedback. – Luke Collins May 23 '18 at 17:42
• Thank you very much for all your optimisations, I have implemented all of them, together with some others (such as replacing the cost lookup table with an numpy matrix). Unfortunately however, I still haven't managed to decrease runtime to under a minute (you said it took you 28 seconds). Is there something I'm missing here? This is the updated code with all changes: ideone.com/oedbmG. – Luke Collins May 24 '18 at 15:03
• You're probably missing the part where that 28 seconds was for a 14-city tour, not a 52-city one. – Peter Taylor May 24 '18 at 15:13
• Oh right, sorry about that. Is it normal for such algorithms to take this long though? Like I said, my friend managed to compile in around 5 seconds in Java. – Luke Collins May 24 '18 at 15:15

An extended comment rather than an answer. You can optimize the cost function as follows:

def cost(self):
if self.__cost is None:
pts = [self.g.vertices[i] for i in self.vertices]
pts.append(self.g.vertices[0])
self.__cost = sum(map(g.d, pts, pts[1:])) + g.d(pts[0], pts[-1])
return self.__cost


First, you avoid looking up each point in g.vertices twice. More importantly, each of [0] + self.vertices and self.vertices + [0] creates a new list, necessarily making a copy of self.vertices. This shaves of about 15% of runtime on my machine.

• Also, this introduces a bug. The original code uses [0] + self.vertices and self.vertices + [0] because 0 is not in self.vertices, but does need to be taken into account in the tour. – Peter Taylor May 23 '18 at 15:03
• @PeterTaylor dist I think is referring to graph.d, the metric function. And as you rightly pointed out, the [0] is excluded for uniqueness of representation (otherwise tours can have multiple representations: each tour is assumed to start from zero) – Luke Collins May 23 '18 at 15:14
• @LukeCollins you are right, I was referring to graph.d. Also, since computing Euclidean distance is pretty cheap, memoizing the distances does save much. It ended up being only a tiny bit faster than recomputing the distance each time. – hilberts_drinking_problem May 23 '18 at 15:53
• @PeterTaylor absolutely right, I missed the point about 0 being present in every tour. – hilberts_drinking_problem May 23 '18 at 15:57
• Given that the point is to save list creation, IMO a better fix is the one I've given in my answer; although reducing from two list creations to one is an improvement. @LukeCollins, for true uniqueness of representation you should enforce t.vertices[0] <= t.vertices[-1] for a tour, since following the tour backwards doesn't fundamentally change it. – Peter Taylor May 23 '18 at 16:01