Simulated Annealing in Python in the early aughties

tsalesp.py

#!/usr/bin/env python

import random
from math import *

# This is an attempt at solving the travelling salesperson problem
# with simulated annealing.  I tried this once in matlab and it sucked.

# a map is a list of tuples, one tuple for each city.
# each tuple represents the (x,y) location of a city on the map.
def create_map(mapsize):
    return [(random.uniform(0,1), random.uniform(0, 1)) for i in range(mapsize)]

# a path is a list of indexes into the list of cities (the map)
# all paths are closed loops and visit each city only once.
def create_random_path(mapsize):
    path = range(mapsize)
    random.shuffle(path)
    return path

# paths are mutated by swapping around city indexes.
# this returns a mutated version of the input path.
def mutate_path(path):
    newpath = path[:]
    index = random.randint(0, len(path)-1)
    newpath[index-1] = path[index]
    newpath[index] = path[index-1] # yay for wraparound arrays  
    return newpath

# mutate a path by randomly selecting a sequence and reversing it
def mutate_path2(path): 
    index1 = random.randint(0, len(path)-2)
    index2 = random.randint(index1+1, len(path)-1)
    segment = path[index1:index2]
    segment.reverse()
    return path[:index1] + segment + path[index2:]

# there's probably a lib function for this,
# but it's faster to write one than to look it up.
def average(list):
    sum = 0.0
    for x in list:
        sum += x
    return sum / len(list)

# euclidian goodness
def distance(ptA, ptB):
    return sqrt((ptB[0] - ptA[0])**2 + (ptB[1] - ptA[1])**2)

# evaluates the length of the given path in the context of the map
def path_length(themap, path):
    totaldist = 0
    for i in range(len(path)):
        totaldist += distance(themap[path[i]], themap[path[i-1]])
    return totaldist

# probability of an energy change based on the Boltzmann distribution
def prob_to_keep(delta_e, temp):
    k = 1.0 # boltzmann's constant in *this* universe
    return exp(-delta_e/(k*temp))

# this is where the cool stuff happens.
# use simulated annealing to find the minimal path length.
def anneal(themap, initialtemp, iterations):
    temp = initialtemp
    tempfactor = 0.95
    path = create_random_path(len(themap))
    consec_successes = 0

    iter = 0  # not using for loop for efficiency reasons
    #while iter < iterations:
    while temp > 0.001:
        # perturb our path a wee bit
        testpath = mutate_path2(path)
        # find the change in energy
        delta_e = path_length(themap, testpath) - path_length(themap, path)

        if delta_e <= 0:
            # ah, this was a step in the Right Direction
            path = testpath
            consec_successes += 1
        else:  # okay, we regressed.  but possibly keep it anyways.
            if random.random() < prob_to_keep(delta_e, temp):
                #print "randomly keeping a crappy move. temp = %f; delta_e = %f." % (temp, delta_e)
                path = testpath
            consec_successes = 0

        # is it time to cool down?
        if consec_successes >= 10 or iter % (100*len(themap)) == 0:
            temp = tempfactor * temp
            consec_successes = 0

        # debug output
        if iter % 1000 == 0:
            print 'DEBUG: iteration %d  temp = %f  bestlen = %f' % (iter, temp, path_length(themap, path))

        iter += 1

    return path


if __name__ == '__main__':
    themap = create_map(10)
    print 'Annealing solutions for map of size %d.' % len(themap)

    bestpath = anneal(themap, 100, 10000)
    print 'The best path had length %f.' % path_length(themap, bestpath)
    avg_rand = average([path_length(themap, create_random_path(len(bestpath))) for i in range(100)])
    print 'Avg rand path has length %f.' % avg_rand
    print 'map used:'
    print themap
    print 'path taken:'
    print bestpath