tsp stuff

WineTourEnv/lib/python2.7/site-packages/django_root/wines/greedy_tsp.py

@@ -0,0 +1,213 @@
+################################################################################
+# A simple algorithm for solving the Travelling Salesman Problem
+# Finds a suboptimal solution
+################################################################################
+try:
+    import psyco
+    psyco.full()
+except:
+    pass
+
+from collections import deque
+
+#logfile = open("greedy.log","w")
+def optimize_solution( distances, connections ):
+    """Tries to optimize solution, found by the greedy algorithm"""
+    N = len(connections)
+    path = restore_path( connections )
+    def ds(i,j): #distance between ith and jth points of path
+        return distances[path[i]][path[j]]
+    d_total = 0.0
+    optimizations = 0
+    for a in xrange(N-1):
+        b = a+1
+        for c in xrange( b+2, N-1):
+            d = c+1
+            delta_d = ds(a,b)+ds(c,d) -( ds(a,c)+ds(b,d))
+            if delta_d > 0:
+                d_total += delta_d
+                optimizations += 1
+                connections[path[a]].remove(path[b])
+                connections[path[a]].append(path[c])
+                connections[path[b]].remove(path[a])
+                connections[path[b]].append(path[d])
+
+                connections[path[c]].remove(path[d])
+                connections[path[c]].append(path[a])
+                connections[path[d]].remove(path[c])
+                connections[path[d]].append(path[b])
+                path[:] = restore_path( connections )
+    #print "Done %d optimizations, total reduction %g"%(optimizations, d_total)
+
+def restore_path( connections ):
+    """Takes array of connections and returns a path.
+    connections is array of lists with 1 or 2 elements.
+    These elements are indices of teh vertices, connected to this vertex
+    """
+    #there are 2 nodes with valency 1 - start and end. Get them.
+    start, end = [ idx for idx, conn in enumerate(connections)
+                   if len(conn)==1 ]
+    path = [start]
+    prev_point = None
+    cur_point = start
+    while True:
+        next_points = [pnt for pnt in connections[cur_point] 
+                       if pnt != prev_point ]
+        if not next_points: break
+        next_point = next_points[0]
+        path.append(next_point)
+        prev_point, cur_point = cur_point, next_point
+    return path
+
+def solve_tsp( distances, optim_steps=3 ):
+    """Given a distance matrix, finds a solution for the TSP problem.
+    Returns list of vertex indices"""
+    N = len(distances)
+    for row in distances:
+        if len(row) != N: raise ValueError, "Matrix is not square"
+
+    #State of the TSP solver algorithm.
+    node_valency = [2] * N #Initially, each node has 2 sticky ends
+    #for each node, stores 1 or 2 connected nodes
+    connections = [[] for i in xrange(N)] 
+
+    def join_segments():
+        #segments of nodes. Initially, each segment contains only 1 node
+        segments = [ [i] for i in xrange(N) ]
+
+        def pairs_by_dist():
+            #print "#### generating vertex pairs"
+            pairs = [None] * (N*(N-1)/2)
+            idx = 0
+            for i in xrange(N):
+                i_n = i*N    
+                dist_i = distances[i]
+                for j in xrange(i+1,N):
+                    pairs[idx] = (dist_i[j],i_n + j) #for the economy of memory, store I and J packed.
+                    idx += 1
+            #print "#### sorting pairs by distance"
+            pairs.sort()
+            return pairs
+
+        def nearest_pairs( sorted_pairs ):
+            for d, i_j in sorted_pairs:
+                i = i_j / N
+                if not node_valency[i] : continue
+                j = i_j % N
+                if not node_valency[j] or (segments[i] is segments[j]): 
+                    continue
+                yield i, j
+
+        pairs_gen = iter( nearest_pairs(pairs_by_dist()) )
+        #print "#### joining segments"
+        for itr in xrange(N-1):
+            i,j = pairs_gen.next()
+            node_valency[i] -= 1
+            node_valency[j] -= 1
+            connections[i].append(j)
+            connections[j].append(i)
+            #join the segments
+            seg_i = segments[i]
+            seg_j = segments[j]
+            if len(seg_j) > len(seg_i):
+                seg_i, seg_j = seg_j, seg_i
+                i, j = j, i
+            for node_idx in seg_j:
+                segments[node_idx] = seg_i
+            seg_i.extend(seg_j)
+
+    join_segments()
+
+    for passn in range(optim_steps):
+        #print "Optimization pass", passn
+        optimize_solution( distances, connections )
+
+    path = restore_path( connections )
+    assert( len(path) == N )
+    return path
+
+def solve_tsp_numpy(node_list, distance_matrix, optim_steps=3 ):
+    """Given a distance matrix, finds a solution for the TSP problem.
+    Returns list of vertex indices.
+    Version that uses Numpy - consumes less memory and works faster."""
+    import numpy
+
+    node_list = sorted(node_list)
+    #logfile.write(str(node_list)+"\n")
+
+    distances = [[[y[z] for z in node_list] for y in distance_matrix][x] for x in node_list]
+
+    N = len(distances)
+    for row in distances:
+        if len(row) != N: raise ValueError, "Matrix is not square"
+
+    #State of the TSP solver algorithm.
+    node_valency = numpy.zeros( N, dtype = 'i1' )
+    node_valency += 2 #Initially, each node has 2 sticky ends
+
+    #for each node, stores 1 or 2 connected nodes
+    connections = [[] for i in xrange(N)] 
+
+    def join_segments():
+        #segments of nodes. Initially, each segment contains only 1 node
+        segments = [ [i] for i in xrange(N) ]
+        def pairs_by_dist_np():
+            #print "#### generating vertex pairs, numpy version"
+            pairs = numpy.zeros( (N*(N-1)/2, ), dtype=('f4, i2, i2') )
+
+            idx = 0
+            for i in xrange(N):
+                row_size = N-i-1
+                dist_i = distances[i]
+                pairs[idx:(idx+row_size)] = [ (dist_i[j], i, j)
+                                              for j in xrange(i+1, N) ]
+                idx += row_size
+            #print "#### sorting pairs by distance"
+            pairs.sort(order=["f1"]) #sort by the first field
+            return pairs
+
+        def nearest_pairs_np( sorted_pairs ):
+            for d, i, j in sorted_pairs:
+                if node_valency[i] and \
+                   node_valency[j] and \
+                   not (segments[i] is segments[j]): 
+                    yield i, j
+
+        pairs_gen = iter( nearest_pairs_np(pairs_by_dist_np()) )
+        #print "#### joining segments"
+        for itr in xrange(N-1):
+            i,j = pairs_gen.next()
+            node_valency[i] -= 1
+            node_valency[j] -= 1
+            connections[i].append(j)
+            connections[j].append(i)
+            #join the segments
+            seg_i = segments[i]
+            seg_j = segments[j]
+            if len(seg_j) > len(seg_i):
+                seg_i, seg_j = seg_j, seg_i
+                i, j = j, i
+            for node_idx in seg_j:
+                segments[node_idx] = seg_i
+            seg_i.extend(seg_j)
+
+    join_segments()
+
+    for passn in range(optim_steps):
+        #print "Optimization pass", passn
+        optimize_solution( distances, connections )
+
+    path = restore_path( connections )
+    distance = 0
+    prev = path[0]
+    for x in path[1:]:
+        distance += distances[prev][x]
+        prev = x
+    distance += distances[path[len(path)-1]][path[0]]
+    assert( len(path) == N )
+    path = [node_list[x] for x in path]
+    path2 = deque(path)
+    path2.rotate(len(path)-path.index(0))
+    path = list(path2)
+    #logfile.write(str(path) + " " + str(distance) + "\n")
+    return [path, distance]

WineTourEnv/lib/python2.7/site-packages/django_root/wines/route_planner.py

@@ -3,18 +3,20 @@
 import os
 import urllib2
 from scipy.sparse.csgraph import dijkstra
-
-from models import Wine, WineProducer
+#from models import Wine, WineProducer
+from operator import itemgetter
+from greedy_tsp import solve_tsp_numpy
 
 PROJECT_ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
 METRIC = {
     'km': 1000,
     'm': 1
 }
 
-def get_coordinates(id):
-    """Given the `id` of a WineProducer, return the latitude and 
-    longitude of the WineProducer as a tuple"""
+
+"""def get_coordinates(id):
+    """""""Given the `id` of a WineProducer, return the latitude and 
+    longitude of the WineProducer as a tuple""""""
     response = urllib2.urlopen('http://maps.googleapis.com/maps' + \
 		'/api/geocode/json?address=%s&sensor=false' \
 		% WineProducer.objects.get(id=id).address.replace(' ', '+'))
@@ -25,9 +27,9 @@ def get_coordinates(id):
 	return 
 
 def get_distance_matrix(ids):
-    """Given an array of WineProducer ids, returns a distance matrix
+    """"""Given an array of WineProducer ids, returns a distance matrix
     of the distances between each point, and a dictionary that maps
-    the rows/cols to the name of the WineProducer it represents"""
+    the rows/cols to the name of the WineProducer it represents""""""
     row_to_name = {}
     coordinate_string = ''
     count = 0
@@ -67,4 +69,136 @@ def get_distance_matrix(ids):
 
 def metric_convert(string):
     distance, unit = string.split(' ')
-    return float(distance) * METRIC[unit]
+    return float(distance) * METRIC[unit]"""
+
+def get_shortest_tours(dist_matrix, constraint_matrix, num_results, min_wineries, max_wineries):
+    """Given a distance matrix where the first row/column corresponds to the source, a constraint matrix for the wineries
+	the number of results desired, the minimum number of wineries in a tour, and the maximum number of wineries in a tour,
+	give a list of the shortest tours ordered by increasing distance"""
+
+    num_wineries = len(constraint_matrix)
+    num_constraints = len(constraint_matrix[0])
+    tour_distance = []
+    max_min_distance = 0
+    tolerance = 1
+    tolerance_percent = 0.02
+    non_improvement_to_count_improvement = {}
+    level_factor = 0
+    def add(a, coeff, b):
+	return a + coeff*b
+
+    log = False
+    if log: logfile = open("tour.log","w")
+
+
+    def calculate_tour_distance(possible_tour):
+	tour = [0] + [x+1 for x in possible_tour] # Because the 0 index in the distance matrix is for the source and incoming indices are one too low
+	return solve_tsp_numpy(tour, dist_matrix)
+
+    def made_improvement(possible_tour, indent):
+        global tour_distance, max_min_distance
+        calculated_tour = calculate_tour_distance(possible_tour)
+        len_tour_dist = len(tour_distance)
+        if log: logfile.write(indent + "Distance: " + str(calculated_tour[1]) + "\n")
+        if len_tour_dist == 0:
+            tour_distance.append(calculated_tour)
+            max_min_distance = tour_distance[0][1]
+            return max_min_distance
+        elif len_tour_dist < num_results:
+            old_high = tour_distance[len_tour_dist-1][1]
+            tour_distance.append(calculated_tour)
+            tour_distance = sorted(tour_distance, key=itemgetter(1))
+            #change = tour_distance[len(tour_distance)-1][1] - max_min_distance
+            max_min_distance = tour_distance[len_tour_dist][1]
+            return max_min_distance-old_high
+        elif calculated_tour[1] < max_min_distance:
+            old_high = tour_distance[len_tour_dist-1][1]
+            tour_distance[num_results-1] = calculated_tour
+            tour_distance = sorted(tour_distance, key=itemgetter(1))
+            max_min_distance = tour_distance[len_tour_dist-1][1]
+            return max_min_distance-old_high
+	return False
+
+    def vec_or(vector1, vector2):
+        ret = []
+        for x in range(len(vector1)):
+            ret.append(vector1[x] or vector2[x])
+        return ret
+
+    def satisfies_new_constraints(possible_tour):
+        len_tour = len(possible_tour)
+        new_node = possible_tour[len_tour-1]
+        already_satisfied = [0] * num_constraints
+        for x in range(len_tour-2):
+            already_satisfied = vec_or(already_satisfied, constraint_matrix[possible_tour[x]])
+        return already_satisfied != vec_or(already_satisfied, constraint_matrix[possible_tour[len_tour-1]])
+
+    def meets_constraints(possible_tour):
+	"""Checks whether a given possible tour satisfies all of the constraints"""
+	constraint_vector = [0]*num_constraints
+	for x in possible_tour:
+	    constraint_vector = vec_or(constraint_vector, constraint_matrix[x])
+	for y in constraint_vector:
+	    if not y:
+		return False
+	return True
+
+    def build_tour_list(consecutive_nonimprovements, possible_tour, level, indent, branch_improvement=False, one_deep_from_constraint_satisfaction = False):
+        global tour_distance, tolerance_percent, level_factor
+	len_current_tour = len(possible_tour)
+        if log: logfile.write(indent + "Possible tour: " + str(possible_tour)+"\n" + indent + "Consecutive nonimprovements: " + str(consecutive_nonimprovements) + "\n")
+        if (consecutive_nonimprovements > tolerance + level*(1+level_factor) and len(tour_distance) == num_results) or len_current_tour > max_wineries:
+            return branch_improvement
+        if meets_constraints(possible_tour):
+            if log: logfile.write(indent + "Meets constraints\n")
+            if len_current_tour >= min_wineries:
+                var = made_improvement(possible_tour, indent)
+                if var and (var > 0 or math.fabs(var)/(tour_distance[len(tour_distance)-1][1] + math.fabs(var)) >= tolerance_percent):
+                    if log: logfile.write(indent + "Made improvement: " + str(var) + " " + str(math.fabs(var)/(tour_distance[len(tour_distance)-1][1] + math.fabs(var))) + " " + str(tolerance_percent) + "\n")
+                    if possible_tour[len_current_tour-1] + 1 < num_wineries:
+                        if not one_deep_from_constraint_satisfaction:
+                            build_tour_list(0, possible_tour + [possible_tour[len_current_tour-1] + 1], level+1, indent + "    ", one_deep_from_constraint_satisfaction = True)
+                            build_tour_list(0, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement)
+                        else:
+                            build_tour_list(0, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement, True)
+                    else:
+                        return True
+                else:
+                    if log: logfile.write(indent + "Didn't make improvement\n")
+                    if possible_tour[len_current_tour-1] + 1 < num_wineries:
+                        build_tour_list(consecutive_nonimprovements + 1, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement, one_deep_from_constraint_satisfaction)
+                    else:
+                        return branch_improvement
+            else:
+                if possible_tour[len_current_tour-1] + 1 < num_wineries:
+                    if build_tour_list(0, possible_tour + [possible_tour[len_current_tour-1] + 1], level+1, indent + "    "):
+                        build_tour_list(0, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement)
+                    else:
+                        build_tour_list(consecutive_nonimprovements+1, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement)
+                else:
+                    return branch_improvement
+        elif possible_tour[len_current_tour-1] + 1 < num_wineries:
+            if not satisfies_new_constraints(possible_tour):
+                if log: logfile.write(indent + "New constraint not met\n")
+                build_tour_list(consecutive_nonimprovements, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement)
+            elif len_current_tour < min_wineries or len_current_tour < num_constraints:
+                if build_tour_list(0, possible_tour + [possible_tour[len_current_tour-1]+1], level+1, indent + "    "):
+                    build_tour_list(0, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement)
+                else:
+                    build_tour_list(consecutive_nonimprovements + 1, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement)
+            else:
+                build_tour_list(consecutive_nonimprovements + 1, [possible_tour[x] for x in range(len_current_tour-1)] + [possible_tour[len_current_tour-1] + 1], level, indent + "    ", branch_improvement)
+        else:
+            return branch_improvement
+
+
+    if log: logfile.close()
+
+    build_tour_list(0, [0], 1, "")
+    return_json = {'tours': []}
+    for x in range(len(tour_distance)):
+        element = {}
+        element['route'] = tour_distance[x][0]
+        element['distance'] = tour_distance[x][1]
+        return_json['tours'].append(element)
+    return return_json