#! /usr/bin/env python
#
#   >->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->
#    AntSystemTSP.py, Heiko Stamer <stamer@informatik.uni-leipzig.de>
#
#      Ant System (AS) for the Traveling Salesman Problem (TSP)
#
#       [The Ant System: Optimization by a colony of cooperating agents]
#                 by M. Dorigo, V. Maniezzo, A. Colorni
#   IEEE Transactions on Systems, Man and Cybernetics - Part B, Vol.26-1 1996
#
#        http://stinfwww.informatik.uni-leipzig.de/~mai97ixb
#   >->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->
#
# Copyright (C) 2000-until_the_end_of_the_ants  <Heiko Stamer>
#
#   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 whrandom
import copy
import math
from time import clock

# Saxony City Tour from "J. Apel: Such- und Graphalgorithmen, Uni-Leipzig"
# OPTIMAL TOUR: 442.0
# [L, D, C, Z, G, R, M]

N = 7
CITIES = range(N)

D = [ 
[0, 110, 67, 72, 190, 65, 80],\
[110, 0, 65, 97, 85, 45, 25],\
[67, 65, 0, 30, 150, 60, 55],\
[72, 97, 30, 0, 180, 87, 85],\
[190, 85, 150, 180, 0, 120, 105],\
[65, 45, 60, 87, 120, 0, 20],\
[80, 25, 55, 85, 105, 20, 0],\
]

def dist(i, j):
	return (math.sqrt(((C[i][0] - C[j][0]) ** 2.0) + ((C[i][1] - C[j][1]) ** 2.0)))

def calc_length(tour):
		l = 0.0
		for part in tour:
			l = l + D[part[0]][part[1]]
		return (l)

################################################################################

class Ant:
	"simple ant agent"
		
	def __init__(self, as, start_city):
		self.AS, self.START_CITY = as, start_city
		
	def __repr__(self):
		return " ANT - START_CITY: %s CURRENT_CITY: %s CURRENT_TOUR: %s\n"\
		% (self.START_CITY, self.CURRENT_CITY, self.CURRENT_TOUR)
				
	def moveTo(self, to_city):
		self.ALLOWED.remove(to_city)
		self.CURRENT_TOUR.append((self.CURRENT_CITY, to_city))
		self.CURRENT_CITY = to_city
		
	def choose(self):
		sum = self.AS.sum_transition(self.CURRENT_CITY, self.ALLOWED)
		p, p_j = whrandom.uniform(0, 1), 0.0
		for j in self.ALLOWED:
			p_j = p_j + (self.AS.transition(self.CURRENT_CITY, j) / sum)
			if (p < p_j):
				return j
			
	def search(self):
		self.CURRENT_CITY, self.CURRENT_TOUR = self.START_CITY, []
		self.ALLOWED = range(N)
		self.ALLOWED.remove(self.CURRENT_CITY)
		while (len(self.ALLOWED) > 0):
			self.moveTo(self.choose())
		self.ALLOWED.append(self.START_CITY)
		self.moveTo(self.START_CITY)
		return (self.CURRENT_TOUR[:])

################################################################################
		
class AntSystem:
	"AS (ant-cycle) environment"
	
	def __init__(self, alpha, beta, rho, q):
		self.ALPHA, self.BETA, self.RHO, self.Q = alpha, beta, rho, q
		self.TAU = copy.deepcopy(D)
		self.dTAU = copy.deepcopy(D)
		self.ACTIVE = copy.deepcopy(D)
		self.ANTS = []
		
	def __repr__(self):
		return " ANT SYSTEM - ALPHA: %s BETA: %s RHO: %s Q: %s\n"\
		% (self.ALPHA, self.BETA, self.RHO, self.Q)
		
	def ETA(self, i, j):
		return ( 1.0 / D[i][j])
	
	def sum_sequence(self, sequence):
		return (reduce((lambda x, y: x + y), sequence, 0))
	
	def init_tau_by_value(self, value):
		for i in CITIES:
			for j in CITIES:
				self.TAU[i][j] = value

	def init_tau_by_matrix(self, matrix):
		for i in CITIES:
			for j in CITIES:
				self.TAU[i][j] = matrix[i][j]
				
	def transition(self, i, j):
		return ((self.TAU[i][j] ** self.ALPHA) * (self.ETA(i, j) ** self.BETA))
	
	def sum_transition(self, i, J):
		sum = 0.0
		for j in J:
			sum = sum + self.transition(i, j)
		return (sum)
	
	def clear_update_transitions(self):
		for i in CITIES:
			for j in CITIES:
				self.dTAU[i][j], self.ACTIVE[i][j] = 0.0, 0
				
	def check_active(self, what):
		sum = 0
		for i in CITIES:
			sum = sum + self.sum_sequence(self.ACTIVE[i])
			if (sum > what):
				return (1)
		return (0)

	def add_update_transitions(self, tour, length):
		for part in tour:
			i, j = part[0], part[1]
			# symmetric TSP
			self.dTAU[i][j] = self.dTAU[i][j] + (self.Q / length)
			self.dTAU[j][i] = self.dTAU[j][i] + (self.Q / length)
			self.ACTIVE[i][j], self.ACTIVE[j][i] = 1, 1

	def update_transitions(self):
		for i in CITIES:
			for j in CITIES:
				self.TAU[i][j] = self.RHO * self.TAU[i][j] + self.dTAU[i][j]
	
	# Initalize	uniform
	def init_uniform(self):
		# M = N and uniformly distributed
		for k in CITIES:
			self.ANTS.append(Ant(self, k))
	
	# Initalize	random
	def init_random(self):
		# M = N and randomly distributed
		for k in CITIES:
			self.ANTS.append(Ant(self, whrandom.choice(CITIES)))
			
	def get_tau(self):
		return (self.TAU[:])
				
	def search(self, T):
		best_tour = []
		best_length = 1000000.0
		no_action_runs = 0
		max_no_action_runs = 10

		# do T iterations of ant-cycle algorithm
		for t in range(0, T):
			self.clear_update_transitions()
			for k in self.ANTS:
				tour = k.search()
				tour_length = calc_length(tour)
				if (tour_length < best_length):
					best_tour = tour
					best_length = tour_length
					print "[%s/%s] local best tour (length = %s):\n%s"\
					% (t + 1, T, best_length, best_tour)
				self.add_update_transitions(tour, tour_length)
			self.update_transitions()
			# checking for stagnation behaviour
			if not(self.check_active(2 * N)):
				no_action_runs = no_action_runs + 1
			if (no_action_runs > max_no_action_runs):
				break
				
		print "[%s/%s] global best tour (length = %s):\n%s"\
		% (t + 1, T, best_length, best_tour)
		print "[%s/%s] iterations done" % (t + 1, T)
		
		return (best_tour)

################################################################################		

as = AntSystem(1.0, 5.0, 0.5, 100.0)
as.init_tau_by_value(0.01)
as.init_uniform()
as.search(300)