trjhtr

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)

 # Add to dictionary
 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):
 """Return total edge-cost of tour"""

 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

 def crossover(self, mum, dad):
 """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]
 bob[i] = dad.vertices[i]

 # Fill the remaining position with the other parents' entries
 current_dad_position, current_mum_position = 0, 0

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

 while dad.vertices[current_dad_position] in alice:
 current_dad_position += 1

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

 alice[i] = dad.vertices[current_dad_position]
 bob[i] = mum.vertices[current_mum_position]

 # Return twins
 return Tour(self.g, alice), Tour(self.g, bob)

 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
 mum, dad = self.select_parent(self.k), self.select_parent(self.k)
 alice, bob = self.crossover(mum, dad)

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

 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)

 # Add to dictionary
 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):
 """Return total edge-cost of tour"""

 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]
 bob[i] = dad.vertices[i]

 # Fill the remaining position with the other parents' entries
 current_dad_position, current_mum_position = 0, 0

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

 while dad.vertices[current_dad_position] in alice:
 current_dad_position += 1

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

 alice[i] = dad.vertices[current_dad_position]
 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 dad.vertices[current_dad_position] in skipalice:
 current_dad_position += 1

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

 alice[i] = dad.vertices[current_dad_position]
 bob[i] = mum.vertices[current_mum_position]
 current_dad_position += 1
 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])
 dadxo = set(dad.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]
 bob[start:start] = dad.vertices[start:end+1]

 # Return twins
 return Tour(self.g, alice), Tour(self.g, bob)

 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.

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

def cost(self):
 """Return total edge-cost of tour"""
 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.