#!/usr/bin/env python
"""
A dirty implementation of the travelling sales man problem as given in the
file siman/siman_tsp.c.
Original Copyright:
/* siman/siman_tsp.c
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
"""
import copy
import Numeric
import pygsl.siman as siman
import pygsl.rng as rng
sin = Numeric.sin
cos = Numeric.cos
pi = Numeric.pi
class City:
def __init__(self, name, lat, longitude):
self._name = name
self._lat = lat
self._longitude = longitude
self._number = None
def __str__(self):
t = (self._name, self._lat, self._longitude, self._number)
return "City % 15s @ (%6.2f, %6.2f) Number %2d" % t
def GetData(self):
return self._lat, self._longitude
def GetName(self):
return Name
def SetNumber(self, number):
self._number = int(number)
def GetNumber(self):
assert(self._number != None)
return self._number
# in this table, latitude and longitude are obtained from the US
# Census Bureau, at http://www.census.gov/cgi-bin/gazetteer */
cities =(City("Santa Fe", 35.68, 105.95),
City("Phoenix", 33.54, 112.07),
City("Albuquerque", 35.12, 106.62),
City("Clovis", 34.41, 103.20),
City("Durango", 37.29, 107.87),
City("Dallas", 32.79, 96.77),
City("Tesuque", 35.77, 105.92),
City("Grants", 35.15, 107.84),
City("Los Alamos", 35.89, 106.28),
City("Las Cruces", 32.34, 106.76),
City("Cortez", 37.35, 108.58),
City("Gallup", 35.52, 108.74))
n_cities = len(cities)
distance_matrix = Numeric.zeros((n_cities, n_cities), Numeric.Float)
cities_vec = Numeric.array(map(lambda x: x.GetData(), cities))
for i in range(len(cities)):
cities[i].SetNumber(i)
earth_radius = 6375.000
# distance between two cities
def city_distance(c):
la = c[:,0]
lo = c[:,1]
sla = sin(la*pi/180)
cla = cos(la*pi/180)
slo = sin(lo*pi/180)
clo = sin(lo*pi/180)
x = cla * clo
y = cla * slo
z = sla
d = Numeric.product(x) * Numeric.product(y) * Numeric.product(z)
return Numeric.arccos(d) * earth_radius
class Tsp:
"""
Travelling sales person
Class describing the problem
"""
def __init__(self):
# Just a sequence of numbers describing the route.
self._route = None
# Just to list the additional required objects
#self.n_cities
#self.r_cities
def GetRoute(self):
return self._route
def SetRoute(self, route):
self._route = route
def SetNCities(self, n_cities):
self.n_cities = n_cities
self.r_cities = Numeric.arange(n_cities)
def EFunc(self):
E = 0
route = self._route
for i in self.r_cities:
E += distance_matrix[route[i]][route[(i + 1) % self.n_cities]]
return E
def Metric(self, other):
route1 = self._route
route2 = other.GetRoute()
distance = 0
for i in self.r_cities:
if route1[i] == route2[i]:
countinue
distance += 1
return distance
def Step(self, rng, step_size):
x1 = rng.get() % (self.n_cities-1) + 1
while 1:
x2 = rng.get() % (self.n_cities-1) + 1
if x1 != x2:
break
x1 = int(x1)
x2 = int(x2)
self._route[x1], self._route[x2] = self._route[x2], self._route[x1]
def Print(self):
for i in self.r_cities:
print " %d " % self._route[i],
def Clone(self):
clone = copy.copy(self)
clone.SetRoute(copy.copy(self._route))
return clone
def prepare_distance_matrix():
for i in range(n_cities):
for j in range(n_cities):
if i == j:
dist = 0
else:
cities = Numeric.take(cities_vec, (i,j))
dist = city_distance(cities)
distance_matrix[i][j] = dist
def print_distance_matrix():
for i in range(n_cities):
for j in range(n_cities):
print "% 8.2f" % distance_matrix[i][j],
print
def siman_exp():
prepare_distance_matrix()
# set up a trivial initial route */
print("# initial order of cities:\n")
for city in cities:
print city
#print " # distance matrix is:"
#print_distance_matrix();
mytsp = Tsp()
mytsp.SetRoute(Numeric.arange(n_cities))
mytsp.SetNCities(n_cities)
result = siman.solve(rng.rng(), mytsp)
route = result.GetRoute()
print("# final order of cities:\n")
for i in route:
print cities[i]
class ExhaustiveSearch:
"""
James Theiler's recursive algorithm for generating all routes
[only works for 12] search the entire space for solutions
"""
def __init__(self):
self.best_E = 1.0e100
self.second_E = 1.0e100
self.third_E = 1.0e100
self.best_route = Numeric.arange(n_cities)
self.second_route = Numeric.arange(n_cities)
self.third_route = Numeric.arange(n_cities)
self.tsp = Tsp()
self.tsp.SetNCities(n_cities)
def do_all_perms(self, route, n):
if n == (n_cities -1):
print "Calculating Energy ...",
self.tsp.SetRoute(route)
self.E = self.tsp.EFunc()
print self.E
if self.E < self.best_E:
self.third_E = self.second_E
self.third_route[:] = self.second_route[:]
self.second_E = self.best_E
self.second_route[:] = self.best_route[:]
self.best_E = self.E
self.best_route[:] = route[:]
elif self.E < self.second_E:
self.third_E = self.second_E
self.third_route[:] = self.second_route[:]
self.second_E = self.E
self.second_route[:] = route[:]
elif self.E < self.third_E:
self.third_E = self.E;
route[:] = self.third_route
else:
new_route = copy.copy(route)
for j in range(n, n_cities):
new_route[j], new_route[n] = new_route[j], new_route[n]
self.do_all_perms(new_route, n+1)
def search(self):
initial_route = Numeric.arange(n_cities)
self.do_all_perms(initial_route, 1)
for i in initial_route:
print cities[i]
print "Best route: "
for i in self.best_route:
print cities[i]
print
print "Second route: "
for i in self.second_route:
print cities[i]
print
print "Third route: "
for i in self.third_route:
print cities[i]
if __name__ == "__main__":
siman_exp()
#ExhaustiveSearch().search()
syntax highlighted by Code2HTML, v. 0.9.1