Sign up to create your own snipts, or login.

brianglass's snipts » Traveling Salesman Problem

posted on Nov 10, 2009 at 8:49 a.m. EST in 
  • import sys
    from django.contrib.gis.db import models
    from django.contrib.localflavor.us import models as localModels
    
    from rmn.settings import DEFAULT_PROJECTION
    
    class Region(models.Model):
    	shape = models.MultiPolygonField( srid=DEFAULT_PROJECTION, null=True, blank=True ) # WGS84 / UTM zone 13N
    	shape_length = models.FloatField( null=True, blank=True )
    	shape_area = models.FloatField( null=True, blank=True )
    	objects = models.GeoManager()
    
    	class Meta:
    		abstract = True
    
    class City(Region):
    	fips = models.IntegerField( null=True, blank=True )
    	name = models.CharField( max_length=30 )
    	url = models.URLField( verify_exists=True, null=True )
    	phone = localModels.PhoneNumberField( null=True )
    
    	def __str__(self): return self.name.title()
    
    	class Meta:
    		ordering = ['name']
    
    class County(Region):
    	fips = models.IntegerField( null=True, blank=True )
    	name = models.CharField( max_length=30 )
    	url = models.URLField( verify_exists=True, null=True )
    	phone = localModels.PhoneNumberField( null=True )
    
    	def __str__(self): return self.name.title()
    
    	class Meta:
    		ordering = ['name']
    
    class Zipcode(Region):
    	code = models.CharField(max_length=5)
    	def __str__(self): return self.code
    
    	class Meta:
    		ordering = ['code']
    
    class Location(models.Model):
    	address1 = models.CharField( max_length=100 )
    	address2 = models.CharField( max_length=100, null=True, blank=True )
    	city = models.ForeignKey( City )
    	state = localModels.USStateField( default='CO' )
    	zip = models.ForeignKey(Zipcode)
    	location = models.PointField( blank=True, srid=DEFAULT_PROJECTION, help_text='The point will be created for you automatically from the address. You may have to click on the map to set a bogus point in order to get this item to save, but it does not matter where the point is.' ) # WGS84 / UTM zone 13N
    	objects = models.GeoManager()
    
    	class Meta:
    		abstract = True
    
    def optimalCircuit(start, destinations):
    	'''
    	Calculate an optimal or near optimal route from the start through each
    	destination point. If the destination list length is less than 15 we do a
    	complete search of the graph. If the list length is larger then we use the greedy
    	algorithm, divide the results in two, and run a full search on the halves.
    
    	Arguments:
    		* start -- a django.contrib.gis.geos.Point object
    		* destinations -- a list of rmn.geo.Location (or a subclass thereof) objects
    
    	'''
    	# precompute distances
    	points = [start] + [ d.location for d in destinations ]
    	distances = [ [0.0] * len(points) for i in range(len(points)) ]
    	for i in range(len(points)):
    		for j in range(i+1, len(points)):
    			distances[i][j] = distances[j][i] = points[i].distance(points[j])
    
    	# greedy algorithm
    	def greedySearch(start, destinations):
    		if len(destinations) > 1:
    			d = [ (distances[start][d],i) for i,d in enumerate(destinations) ]
    			(distance, location) = min( d, key=lambda x:x[0] )
    			# recursion
    			return (destinations[location],) + greedySearch( destinations[location], destinations[:location] + destinations[location+1:] )
    		else:
    			return destinations
    
    	# full search
    	optMemo = {}
    	def fullSearch(start, destinations, end):
    		memokey = (start, destinations, end)
    
    		if len(destinations) > 1:
    			try:
    				# return memoized answer
    				return optMemo[memokey]
    			except KeyError:
    				routes = []
    				for i,d in enumerate(destinations):
    					distance = distances[start][destinations[i]]
    					newDestinations = destinations[:i] + destinations[i+1:]
    					# recursion
    					subpath, subpathLength = fullSearch( d, newDestinations, end )
    					routes.append( ((d,) + subpath, distance + subpathLength) )
    
    				bestRoute, length = min( routes, key=lambda r:r[1] )
    
    				# memoize
    				optMemo[memokey] = (bestRoute, length)
    				return bestRoute, length
    		else:
    			return destinations, distances[start][destinations[0]] + distances[destinations[0]][end]
    
    	# Choose an algorithm
    	if len(destinations) < 15:
    		bestRoute, length = fullSearch(0, tuple(range(1,len(points))), 0)
    		return [ destinations[i-1] for i in bestRoute ]
    	else:
    		greedyRoute = greedySearch(0, tuple(range(1,len(points))))
    		half = len(greedyRoute) // 2
    		bestRoute1, length1 = fullSearch(0, greedyRoute[:half], greedyRoute[half])
    		bestRoute2, length2 = fullSearch(greedyRoute[half-1], greedyRoute[half:], 0)
    		decentRoute = bestRoute1 + bestRoute2
    		return [ destinations[i-1] for i in decentRoute ]
    
    
    # Make circuit finding go faster using the psyco JIT
    try:
    	import psyco
    	psyco.bind(optimalCircuit)
    except ImportError:
    	pass
    

    copy | embed

0 Comments

Sign up, or login to leave a comment.