# coding: utf-8
import matplotlib.pyplot as plt
import numpy as np
import random

########### PARAMETERES OF THE MAP ##############
# Number of cities:
n = 10
# Range of coordinates of cities:
RangeLow = 0
RangeHigh = 10

########### PARAMETERES OF HOPFIELD NETWORK ##############
a = b = 100
c = 50
d = 9
sigma = 2
alpha = 50

########### FUNCTIONS OF HOPFIELD NETWORK ##############
def generateMap():
	# Cities spread on the map randomly
	cities = np.random.uniform(low=RangeLow, high=RangeHigh, size=(n,2))

	# Calculating Euclidan distance between cities using the builtin complex number features
	z=np.array( [[complex(cities[i,0],cities[i,1]) for i in range(n)]] )
	D = abs(z.T-z)
	
	return [cities, D]

def calculateEnergy(D, V):
	# CONSTRAINT #1: Columns should be orthogonal
	dotProduct = V.T.dot(V) 	# Get dot product of each column vectors against each column vector
	dotProduct.flat[::n+1] = 0 # Set main diagonal to zero
	result = a * np.sum( dotProduct )
	
	# CONSTRAINT #2: Rows should be orthogonal
	dotProduct = V.dot(V.T)		# Get dot product of each row vectors against each row vector
	dotProduct.flat[::n+1] = 0 # Set main diagonal to zero
	result += b * np.sum( dotProduct )
	
	# CONSTRAINT #3: Weight of full matrix should be close to n^2
	result += c * ( round(np.sum(V)) - (n + sigma) ) ** 2
	
	# CONSTRAINT #4: Minimize length of path
	result += d * np.sum([V[x,i] * D[x,y] * (V[y,(i+1)%n] + V[y,(i-1)%n])
							for x in range(n) for y in range(n) for i in range(n) if x != y ])
	return result
	
# Activation function using tanh function
def activationFunction(value):
	result = (1 + np.tanh(alpha * value)) / 2
	if result == 0.5: result = 1 # To prevent elements of V matrix to get stuck in the state of 0.5
	return result

# Alternative activation function using sign function
# def activationFunction(value):
	# return 0 if value < 0 else 1
	
def calculateOutputOfNeuron(D, V, x, i):
	return activationFunction	(
		- a * ( np.sum(V[x,:]) - V[x,i] ) 								# CONSTRAINT #1: Columns should be orthogonal
		- b * ( np.sum(V[:,i]) - V[x,i] ) 								# CONSTRAINT #2: Rows should be orthogonal
		- c * ( np.sum(V[:,:]) - (n + sigma) ) 						# CONSTRAINT #3: Weight of full matrix should be close to n^2
		- d * np.sum( D[x,:].dot(V[:,(i+1)%n] + V[:,(i-1)%n]) ) 	# CONSTRAINT #4: Minimize length of path
										)

def optimize(cities, D):
	# INITIALIZATION
	sameEnergyCounter = 0
	cycleCounter = 0
	random.seed()
	# Rows of V represent cities, columns represent the step number
	V = np.random.random((n, n)) 	# Starting form a random state
	# V = np.zeros((n, n))  		# I could start from non-random states as well...
	prevEnergy = calculateEnergy(D, V)

	# OPTIMIZATION
	while cycleCounter < CYCLE_LIMIT and sameEnergyCounter < SAME_LIMIT:
		cycleCounter += 1
		for j in range(5 * n**2):
			x = random.randrange(n)
			i = random.randrange(n)
			V[x,i] = calculateOutputOfNeuron(D, V, x, i)
		energy = calculateEnergy(D, V)
		# For debugging:
		#if cycleCounter % 50 == 0: 
			# print "{:3d}. cycle: {:6.2f}".format(cycleCounter, energy)
			# printTable(V)
		if abs(prevEnergy - energy) < epsilon:
			sameEnergyCounter += 1
		else:
			sameEnergyCounter = 0
		prevEnergy = energy
		#alpha += 1 # It might be a good idea to increase alpha during optimization, but it need some research on how fast it should be increased
		
	return [sameEnergyCounter == SAME_LIMIT, energy, cycleCounter, V]
	
########### HELPING FUNCTIONS ##############
	
# Print a matrix in a nice way: without brackets and rounded up to 1 or down to 0
def printTable(M):
	from prettytable import PrettyTable
	p = PrettyTable()
	for row in M:
		p.add_row([int(round(val,0)) for val in row])
	print p.get_string(header=False, border=False)

# Check if array fulfills constraints
def validate(V):
	# CONSTRAINT #1: Columns should be orthogonal
	dotProduct = V.T.dot(V)
	dotProduct.flat[::n+1] = 0 # Set main diagonal to zero
	result = np.sum( dotProduct )
	
	# CONSTRAINT #2: Rows should be orthogonal
	dotProduct = V.dot(V.T)
	dotProduct.flat[::n+1] = 0 # Set main diagonal to zero
	result += np.sum( dotProduct )
	
	# CONSTRAINT #3: Weight of full matrix should be close to n^2
	result += ( round(np.sum(V)) - n ) ** 2
	
	return result == 0
	
def printDebug(success, energy, cycleCounter, V):
	print "Finishing energy level:", energy
	if success:
		print "Unsuccessful optimalization: cycle limit has reached."
	else:
		print "Successful optimalization after", cycleCounter, "cycles"
		
	print printTable(V)
	print "Does the result fulfills the constraints?", validate(V)
	
# Saves and display the plot of cities and path to the file fn
def plot(cities, C, fn):
	f = plt.figure(1)
	f.clear()
	plt.plot(cities[:,0], cities[:,1], 'o')
	orderedX = cities[:,0].dot(V)
	orderedY = cities[:,1].dot(V)
	plt.plot(np.append(orderedX, orderedX[0]), np.append(orderedY, orderedY[0]))
	f.canvas.draw()
	plt.show(block=False)
	f.savefig(fn)
	return f

########### PARAMETERES OF TESTING ALGORITHM ##############
epsilon = 1 				# If the difference between two numbers are smaller then epsilon, then I consider then equal
CYCLE_LIMIT = 200			# Upper limit of cycles
SAME_LIMIT = 5				# If the energy value holds SAME_LIMIT number of cycles then I assume that I reached local minima
NUMBER_OF_MAPS = 10		# Test should be done on how many maps
NUMBER_OF_TRIALS = 10	# How many times it should try to find the shortest path on a given map

########### TESTING ##############
for mapIndex in range(NUMBER_OF_MAPS):
	print 'Map #' + str(mapIndex)
	[cities, D] = generateMap()
	energyMin = minPlace = -1
	energySum = cycleSum = succesfulTrial = 0
	for trial in range(NUMBER_OF_TRIALS):
		[success, energy, cycleCounter, V] = optimize(cities, D)
		if success and validate(V):
			succesfulTrial += 1
			energySum += energy
			cycleSum += cycleCounter
			f = plot(cities, V, 'result_map' + str(mapIndex) + '_trial' + str(trial))
			if energyMin == -1 or energy < energyMin:
				minPlace = trial
				energyMin = energy
				[energyMin, cycleCounterMin] = [energy, cycleCounter]
	if succesfulTrial > 0:
		print '\tBest result at trial #{:1d} > Length: {:6.3f}, cycle count: {:2d}'.format(minPlace, energyMin, cycleCounterMin)
		print '\tAverage length: {:6.3f}, average cycle count: {:6.3f}'.format(energySum/NUMBER_OF_TRIALS, float(cycleSum)/NUMBER_OF_TRIALS)
		print '\tNumber of succesful trials:', succesfulTrial, '/', NUMBER_OF_TRIALS
	else:
		print '\tNo successful trials are made.'