Tackling the TSP with a genetic algorithm

Question

This is my first go at genetic algorithms, using PyGame for the visual representation. Is there a way to improve the speed of the code? Click to place a point, press ↑ to start looking for the shortest route, or press c to clear all.

import pygame, time, random
#Size of next generation, change number at the end if needed
resize_rate = 1.005
#Controls crossover appearence
cross_rate = 0.7
#Controls number of chromosomes created at the start, exponential
exponent = 4
black = (0, 0, 0)
white = (255, 255, 255)
red = (255, 0, 0)
genes = []
def draw_background(screen):
 screen.fill(white)
def draw_connection(screen, pointlist):
 pygame.draw.lines(screen, black, True, pointlist, 1)
def draw_best_connection(screen, pointlist):
 pygame.draw.lines(screen, red, True, pointlist, 2)
def create_random_chromosomes(genes):
 temp_chromosomes = []
 chromosomes = []
 for i in range(len(genes)**exponent):
 genes_copy = genes[:]
 random.shuffle(genes_copy)
 temp_chromosomes.append(genes_copy)
 for element in temp_chromosomes:
 if element not in chromosomes:
 chromosomes.append(element)
 print(str(len(chromosomes)) + " random chromosomes created")
 return chromosomes
def calculate_distance(chromosome):
 distances = []
 total_distance = 0
 for i in range(len(chromosome)):
 if i < (len(chromosome)-1):
 distance_A_B = (((chromosome[i+1][0]-chromosome[i][0])**2)+(chromosome[i+1][1]-chromosome[i][1])**2)**0.5
 distances.append(distance_A_B)
 else:
 distance_A_B = (((chromosome[i][0]-chromosome[0][0])**2)+(chromosome[i][1]-chromosome[0][1])**2)**0.5
 distances.append(distance_A_B)
 for i in range(len(distances)):
 total_distance += distances[i]
 return total_distance
def calculate_fitness(chromosomes):
 fitness_scores = []
 for i in range(len(chromosomes)):
 fitness = 0
 fitness = 1/calculate_distance(chromosomes[i])
 fitness_scores.append(fitness)
 return fitness_scores
def build_chromosome_fitness_dict(chromosomes, fitness_scores):
 chrom_fit_dict = {}
 for i in range(len(chromosomes)):
 chrom_fit_dict[fitness_scores[i]] = chromosomes[i]
 return chrom_fit_dict
def roulette_selection(chromosomes, fitnesses, new_size):
 sum_fitness = sum(fitnesses)
 rel_fitness = [fitness/sum_fitness for fitness in fitnesses]
 #Generate probability intervals for each chromosome
 probabilities = [sum(rel_fitness[:i+1]) for i in range(len(rel_fitness))]
 #Select chromosomes
 selected_chromosomes = []
 for n in range(new_size):
 r = random.random()
 for (i, chromosome) in enumerate(chromosomes):
 if r <= probabilities[i]:
 selected_chromosomes.append(chromosome)
 break
 return selected_chromosomes
def find_closest(cities):
 dist_a = (((cities[0][0]-cities[1][0])**2)+(cities[0][1]-cities[1][1])**2)**0.5
 dist_b = (((cities[0][0]-cities[2][0])**2)+(cities[0][1]-cities[2][1])**2)**0.5
 if dist_a < dist_b:
 return cities[1]
 else:
 return cities[2]
#Greedy crossover selects the first city of one parent,
#compares the cities leaving that city in both parents,
#and chooses the closer one to extend the tour.
#If one city has already appeared in the tour, we choose the other city.
#If both cities have already appeared, we randomly select a non-selected city. 
def crossover(chromosomes, rate):
 new_chromosomes = []
 #len -1 because index out of range error might occur otherwise
 for i in range(len(chromosomes)-1):
 #Crossover T/F
 if random.random() < rate:
 #Take 2 chromosomes
 a = chromosomes[i]
 b = chromosomes[i+1]
 #Randomly select crossover location
 cross_location = random.randint(1, len(chromosomes[i]))
 #Create parts
 c = a[:cross_location]
 for i in range((len(a)-len(c))):
 #City to start from
 start_pos = c[-1]
 #2 most plausible targets
 target_0 = a[cross_location+i]
 target_1 = b[cross_location+i]
 #Determine if targets already appear in c and add correct city
 if target_0 in c:
 if target_1 in c:
 not_yet_selected = []
 for city in genes:
 if city not in c:
 not_yet_selected.append(city)
 c.append(not_yet_selected[random.randint(0,(len(not_yet_selected)-1))])
 else:
 c.append(target_1)
 elif target_1 in c:
 c.append(target_0)
 else:
 cities = [start_pos, target_0, target_1]
 closest = find_closest(cities)
 c.append(closest)
 new_chromosomes.append(c)
 else:
 new_chromosomes.append(chromosomes[i])
 return new_chromosomes
pygame.init()
screen = pygame.display.set_mode((640,340))
clock = pygame.time.Clock()
done = False
first = True
while done == False:
 for event in pygame.event.get():
 if event.type == pygame.QUIT:
 done = True
 if event.type == pygame.KEYDOWN:
 if event.key == pygame.K_UP:
 if genes != []:
 #Start Genetic Algorithm
 print('Creating initial solutions')
 chromosomes = create_random_chromosomes(genes)
 #Draw most possible solutions
 for i in range(len(chromosomes)):
 draw_connection(screen, chromosomes[i])
 #Main Logic Loop
 print('Starting calculations')
 while len(chromosomes) > 10:
 print(str(len(chromosomes)) + ' chromosomes left')
 new_size = int(len(chromosomes)/resize_rate)
 #Calculate fitness scores
 fitness_scores = calculate_fitness(chromosomes)
 #Build dictionary to preserve order
 chrom_fit_dict = build_chromosome_fitness_dict(chromosomes, fitness_scores)
 #Roulette selection - select chromosomes for next generation
 selected_chromosomes = roulette_selection(chromosomes, fitness_scores, new_size)
 #Adopt newly selected chromosomes as population
 chromosomes = selected_chromosomes
 #Apply crossover mutation
 mutated_chromosomes = crossover(chromosomes, cross_rate)
 #Adopt mutated chromosomes as population
 chromosomes = mutated_chromosomes
 best_chrom = fitness_scores.index(max(chrom_fit_dict))
 draw_best_connection(screen, chromosomes[best_chrom])
 if event.key == pygame.K_c:
 draw_background(screen)
 genes = []
 if first == True:
 draw_background(screen)
 first = False
 #Get points
 if pygame.mouse.get_pressed() == (1, 0, 0):
 mouse_pos = pygame.mouse.get_pos()
 genes.append(mouse_pos)
 pygame.draw.rect(screen, black, [mouse_pos[0], mouse_pos[1], 2, 2], 1)
 genes = list(set(genes))
 pygame.display.flip()
 clock.tick(30)
pygame.quit()

• cat_baxter: How many vertexes (cities) do you want to process? Do you need a genetic algorithm only, not an exact one?
• Daquicker: It isn't meant for real-world application so it doesn't really matter how many cities/points it can process, largest test I ran was with 30 different points, it executed in about 7-8 minutes and was most of the times I ran it 90% - 100% correct. I'm only interested in the genetic algorithm, this was a test project to learn how these work, so an exact result is not requiered.
• cat_baxter: Ok, I see. If you need the exact TSP solution (may be in order to check you genetic algorithm solution) let me know. It solves 30 points TSP in about 0.5 second.
• Daquicker: oh, to check solutions that would be interesting indeed.

Answer 1

Ok, so here is exact ATSP solver, just in order to check your genetic algorithm solution. It's based on pascal code implementing the branch-and- bound algorithm by Little, Murty, Sweeney, and Karel. It's very very fast (30 points in 0.5 sec!).

'''
Branch and bound Algorithm for Travelling Salesman Problem
based on
http://www.cs.sunysb.edu/~algorith/implement/syslo/distrib/processed/babtsp.p
'''
import time
INF = 100000000
best_cost = 0
def reduce(size, w, row, col, rowred, colred):
 rvalue = 0
 for i in xrange(size):
 temp = INF
 for j in xrange(size):
 temp = min(temp, w[row[i]][col[j]])
 if temp > 0:
 for j in xrange(size):
 if w[row[i]][col[j]] < INF:
 w[row[i]][col[j]] -= temp
 rvalue += temp
 rowred[i] = temp
 for j in xrange(size):
 temp = INF
 for i in xrange(size):
 temp = min(temp, w[row[i]][col[j]])
 if temp > 0:
 for i in xrange(size):
 if w[row[i]][col[j]] < INF:
 w[row[i]][col[j]] -= temp
 rvalue += temp
 colred[j] = temp
 return rvalue
def bestEdge(size, w, row, col):
 mosti = -INF
 ri = 0
 ci = 0
 for i in xrange(size):
 for j in xrange(size):
 if not w[row[i]][col[j]]:
 minrowwelt = INF
 zeroes = 0
 for k in xrange(size):
 if not w[row[i]][col[k]]:
 zeroes += 1
 else:
 minrowwelt = min(minrowwelt, w[row[i]][col[k]])
 if zeroes > 1: minrowwelt = 0
 mincolwelt = INF
 zeroes = 0
 for k in xrange(size):
 if not w[row[k]][col[j]]:
 zeroes += 1
 else:
 mincolwelt = min(mincolwelt, w[row[k]][col[j]])
 if zeroes > 1: mincolwelt = 0
 if minrowwelt + mincolwelt > mosti:
 mosti = minrowwelt + mincolwelt
 ri = i
 ci = j
 return mosti, ri, ci
def explore(n, w, edges, cost, row, col, best, fwdptr, backptr):
 global best_cost
 colred = [0 for _ in xrange(n)]
 rowred = [0 for _ in xrange(n)]
 size = n - edges
 cost += reduce(size, w, row, col, rowred, colred)
 if cost < best_cost:
 if edges == n - 2:
 for i in xrange(n): best[i] = fwdptr[i]
 if w[row[0]][col[0]] >= INF:
 avoid = 0
 else:
 avoid = 1
 best[row[0]] = col[1 - avoid]
 best[row[1]] = col[avoid]
 best_cost = cost
 else:
 mostv, rv, cv = bestEdge(size, w, row, col)
 lowerbound = cost + mostv
 fwdptr[row[rv]] = col[cv]
 backptr[col[cv]] = row[rv]
 last = col[cv]
 while fwdptr[last] != INF: last = fwdptr[last]
 first = row[rv]
 while backptr[first] != INF: first = backptr[first]
 colrowval = w[last][first]
 w[last][first] = INF
 newcol = [INF for _ in xrange(size)]
 newrow = [INF for _ in xrange(size)]
 for i in xrange(rv): newrow[i] = row[i]
 for i in xrange(rv, size - 1): newrow[i] = row[i + 1]
 for i in xrange(cv): newcol[i] = col[i]
 for i in xrange(cv, size - 1): newcol[i] = col[i + 1]
 explore(n, w, edges + 1, cost, newrow, newcol, best, fwdptr, backptr)
 w[last][first] = colrowval
 backptr[col[cv]] = INF
 fwdptr[row[rv]] = INF
 if lowerbound < best_cost:
 w[row[rv]][col[cv]] = INF
 explore(n, w, edges, cost, row, col, best, fwdptr, backptr)
 w[row[rv]][col[cv]] = 0
 for i in xrange(size):
 for j in xrange(size):
 w[row[i]][col[j]] = w[row[i]][col[j]] + rowred[i] + colred[j]
def atsp(w):
 global best_cost
 size = len(w)
 col = [i for i in xrange(size)]
 row = [i for i in xrange(size)]
 best = [0 for _ in xrange(size)]
 route = [0 for _ in xrange(size)]
 fwdptr = [INF for _ in xrange(size)]
 backptr = [INF for _ in xrange(size)]
 best_cost = INF
 explore(size, w, 0, 0, row, col, best, fwdptr, backptr)
 index = 0
 for i in xrange(size):
 route[i] = index
 index = best[index]
 index = []
 cost = 0
 for i in xrange(size):
 if i != size - 1:
 src = route[i]
 dst = route[i + 1]
 else:
 src = route[i]
 dst = 0
 cost += w[src][dst]
 index.append([src, dst])
 return cost, index
def main():
 # adjasted matrix
 m = [
 [INF, 514, 230, 92, 172, 201, 320, 205, 329, 285, 404, 128, 352, 357, 190, 269, 484, 237, 567, 209, 413, 404,
 333, 480, 371, 175, 470, 281, 423, 328],
 [452, INF, 359, 410, 450, 401, 315, 293, 244, 200, 275, 446, 469, 143, 278, 362, 399, 344, 544, 644, 335, 381,
 248, 457, 626, 375, 551, 258, 397, 305],
 [414, 360, INF, 90, 170, 42, 90, 46, 282, 238, 186, 126, 320, 198, 31, 267, 437, 190, 358, 623, 246, 345, 286,
 271, 267, 16, 326, 222, 421, 119],
 [324, 427, 269, INF, 80, 109, 359, 113, 263, 233, 312, 36, 260, 265, 98, 177, 418, 257, 475, 533, 321, 312, 281,
 388, 279, 83, 378, 189, 331, 236],
 [244, 382, 312, 226, INF, 335, 402, 246, 197, 153, 334, 262, 486, 466, 231, 315, 352, 297, 670, 453, 241, 538,
 201, 583, 505, 309, 298, 415, 350, 431],
 [448, 318, 355, 216, 268, INF, 445, 111, 240, 196, 288, 252, 385, 156, 96, 358, 395, 148, 691, 640, 348, 528,
 244, 455, 373, 299, 566, 405, 393, 452],
 [330, 437, 126, 216, 280, 119, INF, 123, 206, 315, 312, 252, 348, 275, 108, 324, 361, 267, 268, 539, 372, 255,
 242, 181, 393, 93, 236, 132, 478, 29],
 [495, 557, 521, 539, 313, 563, 469, INF, 406, 444, 543, 575, 274, 700, 440, 524, 442, 506, 729, 633, 554, 716,
 442, 642, 788, 537, 611, 593, 641, 490],
 [208, 256, 115, 168, 206, 157, 205, 49, INF, 148, 137, 204, 323, 313, 34, 118, 155, 100, 473, 417, 197, 460, 36,
 386, 382, 131, 441, 337, 272, 234],
 [252, 229, 159, 212, 250, 201, 249, 93, 44, INF, 181, 248, 367, 357, 78, 162, 199, 144, 517, 444, 241, 398, 48,
 430, 426, 175, 485, 275, 197, 278],
 [534, 568, 363, 453, 411, 250, 237, 254, 326, 344, INF, 225, 449, 406, 239, 366, 481, 398, 505, 547, 60, 375,
 362, 418, 623, 224, 473, 252, 520, 266],
 [445, 507, 334, 389, 469, 332, 393, 336, 413, 369, 276, INF, 224, 488, 321, 141, 392, 321, 659, 583, 336, 496,
 417, 572, 601, 306, 629, 373, 295, 420],
 [221, 283, 313, 313, 393, 314, 195, 263, 214, 170, 351, 349, INF, 426, 248, 325, 168, 314, 455, 359, 411, 442,
 218, 368, 580, 288, 431, 319, 367, 216],
 [309, 286, 216, 269, 307, 258, 306, 150, 101, 57, 132, 305, 326, INF, 135, 219, 256, 201, 574, 501, 192, 455,
 105, 487, 483, 232, 542, 332, 254, 335],
 [416, 554, 484, 398, 172, 507, 484, 15, 369, 325, 506, 434, 289, 638, INF, 487, 457, 469, 717, 625, 413, 554,
 373, 359, 677, 481, 470, 431, 522, 478],
 [514, 509, 193, 283, 352, 191, 252, 195, 306, 387, 379, 319, 469, 347, 180, INF, 461, 339, 518, 723, 249, 355,
 342, 431, 460, 165, 488, 232, 154, 279],
 [53, 115, 278, 145, 225, 254, 368, 212, 163, 119, 300, 181, 405, 258, 197, 281, INF, 263, 620, 262, 360, 457,
 167, 533, 424, 228, 523, 334, 316, 381],
 [300, 277, 207, 68, 148, 177, 297, 141, 92, 48, 229, 104, 328, 333, 126, 210, 247, INF, 543, 492, 289, 380, 96,
 456, 347, 151, 446, 257, 245, 304],
 [758, 471, 620, 710, 754, 593, 710, 597, 550, 671, 687, 746, 571, 614, 582, 668, 705, 650, INF, 930, 747, 599,
 586, 404, 887, 567, 915, 476, 822, 523],
 [39, 434, 269, 131, 211, 240, 359, 227, 178, 324, 315, 167, 391, 396, 212, 296, 333, 276, 452, INF, 375, 439,
 214, 365, 410, 214, 509, 316, 450, 213],
 [526, 513, 443, 496, 534, 485, 533, 377, 328, 284, 441, 165, 389, 641, 362, 306, 483, 428, 698, 487, INF, 535,
 332, 611, 710, 459, 769, 412, 460, 459],
 [784, 654, 691, 552, 604, 336, 781, 447, 576, 532, 624, 588, 721, 492, 432, 694, 731, 484, 1027, 976, 684, INF,
 580, 791, 709, 635, 902, 741, 729, 788],
 [435, 220, 576, 527, 607, 569, 450, 513, 464, 420, 495, 563, 689, 363, 498, 582, 619, 564, 718, 396, 555, 601,
 INF, 631, 806, 543, 686, 478, 617, 479],
 [354, 402, 216, 306, 350, 189, 306, 193, 146, 294, 283, 342, 167, 345, 178, 264, 301, 246, 358, 526, 343, 195,
 182, INF, 483, 163, 511, 72, 418, 119],
 [311, 288, 218, 79, 159, 188, 308, 152, 103, 59, 240, 115, 113, 344, 137, 221, 258, 11, 554, 472, 300, 391, 107,
 462, INF, 162, 457, 268, 256, 315],
 [431, 344, 186, 242, 187, 26, 276, 30, 266, 222, 314, 278, 304, 182, 15, 384, 421, 174, 544, 640, 374, 531, 270,
 374, 399, INF, 348, 408, 419, 305],
 [477, 454, 384, 437, 475, 426, 474, 318, 269, 225, 300, 473, 494, 168, 303, 387, 424, 369, 742, 669, 360, 623,
 273, 655, 651, 400, INF, 500, 422, 503],
 [282, 330, 144, 234, 280, 186, 234, 123, 74, 222, 211, 270, 366, 342, 108, 192, 229, 174, 286, 491, 271, 123,
 110, 199, 411, 160, 470, INF, 346, 47],
 [360, 399, 39, 129, 209, 81, 98, 85, 152, 277, 225, 165, 359, 237, 70, 270, 307, 229, 364, 569, 95, 201, 188,
 277, 306, 55, 334, 78, INF, 125],
 [301, 433, 97, 187, 267, 139, 187, 143, 177, 325, 283, 223, 319, 295, 128, 295, 332, 277, 239, 510, 343, 226,
 213, 152, 364, 113, 423, 103, 449, INF]
 ]
 start_time = time.time()
 cost, path = atsp(m)
 print "Cost = ", cost
 print "Path = ", path
 print "Time (s)", time.time() - start_time
if __name__ == "__main__":
 main()