Tabu.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
import matplotlib.pyplot as plt
import time

CallsLimit = 5000
dimension = 2

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 = [0.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 frange(start, stop, step):
    i = start
    while i < stop:
        yield i
        i += step

def constraintsPlot():
	_plot = [[0.075], [10]]
	for i in frange(0.1, 10.0, 0.1):
		constr = [[i], [0.75/i]]
		_plot = np.concatenate((_plot, constr), axis = 1)
	for i in frange(5.0, 11.0, 1.0):
		constr = [[15.0-i], [i]]
		_plot = np.concatenate((_plot, constr), axis = 1)
	_plot = np.concatenate((_plot, [[0.075], [10]]), axis = 1)
	return _plot

def PlotSearch(_plot):
	plt.plot(zip(*_plot)[0],zip(*_plot)[1], '.')
	_plotConstr = constraintsPlot()
	plt.plot(_plotConstr[0],_plotConstr[1])
	plt.plot(1.6070290311590574, 0.46672081580533736, '*')
	plt.show()
	plt.clf()

def ListUpdate(List, newPoint):
	l = len(List)-1
	for i in range(l):
		List[l-i] = List[l-i-1]
	List[0] = newPoint

def isAllowed(x):
	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):
	values = []
	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	
		
		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 isVisited(x):
	tol = 1.e-8
	for element in steps:
		dist = np.linalg.norm([a-b for a,b in zip(x,element)])
		if (dist < tol):
			return True

	return False

def MakeStep(x, tabuList, current_value):
	values = CheckSteps(x)

	best_neigh_value = np.min(values)
	best_arg = np.argmin(values)
	
	if (best_neigh_value == 1.0):
		return x_init(), objFunction(x_init())

	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 CheckNeighbors(x):
	neighCount = 0
	radius = 1.5#0.6*math.sqrt(2.0*dimension)
	for pos in steps:
		if (np.linalg.norm([a-b for a,b in zip(x,pos)]) < dx*radius):
			neighCount += 1
		if neighCount > 20:#2.0*math.sqrt(4.0/dimension)*5.0:
			return True
	return False

#plt.figure()

dx_init = 4.0
diffStepsValues = []
while (dx_init < 4.2):

	finalValues = []
	for ITER in range(1):
		
		np.random.seed(1)

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

		steps = []

		MTMSize = int(4.0*(1.0+dimension/5))
		BestValues = [1.0*i for i in range(MTMSize)]
		BestPositions = [[0.0 for d in range(dimension)] for i in range(MTMSize)]
		MTMCount = 0; MTMCountMax = 10

		StepMCount = 0
		StepMCountMax = 25#*(1.0+dimension/10.0)

		Dmin = 3.0*dx
		Dsim = 0.5*dx

		while (TotalCalls < CallsLimit):

			if (TotalCalls < CallsLimit*0.9):
				while (CheckNeighbors(x) ): x = x_init()

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

			if (value_new < np.max(BestValues)): #if improvement was made
				if not (isVisited(x_new)):
					if (value_new < np.min(BestValues)):
						#print 'Improvement made: ', value_new, ' Total calls: ', TotalCalls, 'dx: ', dx#, ' at: ', x_new
						StepMCount = 0
						MTMCount = 0

					BestPositions[np.argmax(BestValues)] = x_new
					BestValues[np.argmax(BestValues)] = value_new
			else:
				MTMCount += 1
				StepMCount += 1

			if (MTMCount >= MTMCountMax):
				for d in range(dimension):
					x_new[d] = 0.0
					for j in range(MTMSize):
						x_new[d] += BestPositions[j][d]/MTMSize
				MTMCount = 0
				#print 'Medium term memory activated'

			if (StepMCount >= StepMCountMax):
				x_new = BestPositions[np.argmin(BestValues)]
				value_new = np.min(BestValues)
				#print 'Going to point with ', np.min(BestValues), 'at ', x_new
				#tabuList = [[0.0 for d in range(dimension)] for i in range(7)] #initially no tabus
				if (TotalCalls < CallsLimit*0.7):
					dx = dx*(1-0.3)
				else: dx = dx * (1-0.5)
				Dmin = 3.0*dx; Dsim = 0.5*dx
				if (dx < 1.e-6): break
				StepMCount = 0
				MTMCount = 0

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

		finalValues.append(np.min(BestValues))
		print np.min(BestValues)


	tmpcount = 0
	for val in finalValues:
		if val < -0.7: tmpcount += 1
	print np.mean(finalValues), np.std(finalValues), dx_init, 'success rate: ', tmpcount*100.0/len(finalValues), '%'
	dx_init += 0.2

	diffStepsValues.append(finalValues)

#plt.plot(finalValues)
#plt.show()

if (dimension == 2):
	PlotSearch(steps)

print diffStepsValues