Tabu_show.py

#10D Improvement made:  -0.737295065591  at:  [3.1419642100697103, 3.083081517990897, 3.037046544480138, 2.983031641774718, 1.4454720306407767, 0.4157516943765218, 0.3054390332826271, 0.25557513716674896, 0.3844567203155933, 0.47479136622732576] 
import math
import random
import numpy as np

CallsLimit = 5000
dimension = 10

def objFunction(x):
	global TotalCalls; TotalCalls += 1
	_cos_summ = 0.0	
	_cos_prod = 1.0
	_denomin  = 0.0
	i = 1.0
	for xi in x:
		_cos_summ += (math.cos(xi))**4.0
		_cos_prod *= (math.cos(xi))**2.0
		_denomin += i*xi**2.0
		i += 1.0

	result = (_cos_summ - 2.0 * _cos_prod)/math.sqrt(_denomin)

	return -result
def ConstraintsCheck(x):
	_constr_prod = 1.0
	_constr_sum  = 0.0
	for xi in x:
		if (xi < 0.0 or xi > 10.0):
			return False
		_constr_prod *= xi
		_constr_sum += xi

	if (_constr_prod <= 0.75):
		return False
	if (_constr_sum >= 7.5*len(x)):
		return False
	return True

def x_init():
	x_init = [10.0 for i in range(dimension)]
	while not ConstraintsCheck(x_init):
		x_init = [10.0*np.random.rand() for i in range(dimension)]
	return x_init

def ListUpdate(List, newPoint): #Updating tabu list. First come, first leave
	l = len(List)-1
	for i in range(l):
		List[l-i] = List[l-i-1]
	List[0] = newPoint

def isAllowed(x):				#We do not want to visit tabu locations
	tol = 1.e-8
	for i in range(len(tabuList)):
		dist = np.linalg.norm([a-b for a,b in zip(x,tabuList[i])])
		if (dist < tol):
			return False
	return True

def CheckSteps(x):						 #Check possible moves
	values = []							 #values near the current x position
	values.append([]); values.append([]) #first column for negative direction x-dx, second for positive x+dx

	for direction in range(dimension):

		step = [0.0 for i in range(dimension)]; step[direction] = dx	#step vector
		
		xNew = [a-b for a,b in zip(x,step)] #making step in negative direction
		if (ConstraintsCheck(xNew) and isAllowed(xNew)):		#if step is in the allowed range
			values[0].append(objFunction(xNew))	
		else:							#else, calculate cost function value
			values[0].append(1.0)

		xNew = [a+b for a,b in zip(x,step)] #making step in positive direction
		if (ConstraintsCheck(xNew) and isAllowed(xNew)):		
			values[1].append(objFunction(xNew))	
		else: 
			values[1].append(1.0)

	return values


def MakeStep(x, tabuList, current_value):		#Making best step
	values = CheckSteps(x)

	best_neigh_value = np.min(values)			#The best possible step
	best_arg = np.argmin(values)	
	
	if (best_neigh_value == 1.0):				#If everything is already tabu or constrained
		return x_init(), objFunction(x_init())  #Make a random step

	step = [0.0 for i in range(dimension)]
	step[best_arg%dimension] = dx
	if (best_arg < dimension): #meaning that best step is in negative direction
		direction = -1.0
	else: direction = 1.0

	xNew = [x[i] + direction*step[i] for i in range(dimension)]

	ListUpdate(tabuList,xNew); x = xNew

	if (best_neigh_value < current_value): #if new value is better than current one, make pattern move
		xNew = [x[i] + direction*step[i] for i in range(dimension)]
		if (ConstraintsCheck(xNew)):
			ListUpdate(tabuList,xNew)
		else:
			xNew = x

	return xNew, best_neigh_value


def distance(x,y):			#Distance between two points
	z = [x[i] - y[i] for i in range(dimension)]
	return np.linalg.norm(z)


def D_Condition(d):					   #If there are no neighbors in d range
	for pos in BestPositions:
		if (distance(pos, x_new) < d): #If there is a point in archive in d range
			return False   				  #Then we return False

	return True

def FindNearest():			#The nearest neighbor for archiving
	distances = []
	for pos in BestPositions:
		distances.append(distance(pos, x_new))
	return np.argmin(distances)  

def InitD(): 				#Dsim and Dmin as functions of dx
	return 1.2*dx, 0.12*dx

dx_init = 2.6
SEED = 1
np.random.seed(SEED)

TotalCalls = 0
x = x_init()
current_value = objFunction(x)
dx = dx_init
tabuList = [[0.0 for d in range(dimension)] for i in range(10)] #initially no tabus

steps = [] #Storage of all visited locations

MTMSize = int(4.0*(1.0+dimension/10))
BestPositions = [x_init() for i in range(MTMSize)]
BestValues = [objFunction(BestPositions[i]) for i in range(MTMSize)]

MTMCount = 0; MTMCountMax = 10
StepMCount = 0; StepMCountMax = 25

Dmin, Dsim = InitD()

while (TotalCalls < CallsLimit):

	x_new, value_new = MakeStep(x, tabuList, current_value)

	if (value_new < np.min(BestValues)): #If the new best new solution has been found
		MTMCount = 0; StepMCount = 0 
	else:
		MTMCount+= 1; StepMCount+= 1

	#Smart Archiving
	if (D_Condition(Dmin)):		#There is no solutions close to x_new
		BestPositions[np.argmax(BestValues)] = x_new
		BestValues[np.argmax(BestValues)] = value_new
	else:						#There is solutions close to x_new
		if (value_new < np.min(BestValues)): #We have found the best new solution and there is a close neighbor
			nearest = FindNearest()
			BestPositions[nearest] = x_new
			BestValues[nearest] = value_new
		else:
			if not (D_Condition(Dsim)):			 # We can improve an already existing solution
				nearest = FindNearest()
				if (value_new < BestValues[nearest]):
					BestPositions[nearest] = x_new
					BestValues[nearest] = value_new

	#Medium term Memory trigger
	if (MTMCount >= MTMCountMax):
		x_new = np.mean(BestPositions, axis = 0)
		MTMCount = 0

	#Step reduction trigger
	if (StepMCount >= StepMCountMax):
		x_new = BestPositions[np.argmin(BestValues)]
		value_new = np.min(BestValues)
		if (TotalCalls < CallsLimit*0.9):
			dx = dx*(1-0.35)
		else: dx = dx * (1-0.5)
		Dmin, Dsim = InitD()
		StepMCount = 0
		MTMCount = 0

	steps.append(x_new)
	current_value = value_new
	x = x_new