SA.py

import math
import random
import numpy as np
import time

CallsLimit = 5000
dimension = 10

TrialsNumber = 199
AcceptsNumber = int(0.6*TrialsNumber)

DStepSize = 1.0

alpha = 0.1
omega = 2.1

TempAdj = 20.0

def objFunction(x): #Objective function value for given 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): #Checks if the given x lies inside the given region
	_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 D_init(): #D matrix initialization
	D = [DStepSize for i in range(dimension)]
	return D
def x_init(): #Initial random x position
	x_init = [10.0*np.random.rand() 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 D_upd(D, step): #Updating D matrix after accepted step. Upper and lower limits applied
	Dnew = [min(max((1.0-alpha)*D[i] + alpha*omega*math.fabs(step[i]), 1.e-4), DStepSize) for i in range(dimension)]
	return Dnew

def step(D, u): #Calculate the step dx = Du
	step = [u[i]*D[i] for i in range(dimension)]
	return step

def frange(start, stop, step):
    i = start
    while i < stop:
        yield i
        i += step

def getInitTemp(): #Evaluating initial temperature - accepting all the values, then calculating the standard deviation
	D = D_init()
	x = x_init()
	Values = [] #List of visited positions values
	for i in range(int(0.01*CallsLimit)):
		u = 2.0*(np.random.rand(dimension, 1) - 0.5) #Performing random step
		x = x + step(D,u) 							 #Updating position
		j = 0
		while not ConstraintsCheck(x):               #If not inside the region, doing new random step
			j+=1
			u = 2.0*(np.random.rand(dimension, 1) - 0.5)
			x = x + step(D,u)
			if (j > 10000):                          #If we are stuck at one point for too long, go to another random direction
				x = 10.0*np.random.rand(dimension, 1)
				j = 0

		Values.append(objFunction(x))
	return np.std(Values)

def MakeStep(x, D):
	xnew = [0.0 for i in range(dimension)]
	while not ConstraintsCheck(xnew):
		u = np.random.normal([0.0 for i in range(dimension)], [TempAdj*math.sqrt(T) for i in range(dimension)], dimension)
		#u = [2.0*(np.random.rand()-0.5) for i in range(dimension)]
		xnew = [x[i] + step(D,u)[i] for i in range(dimension)]
	
	return step(D,u)


def T_upd(T, sigma): #Updating temperature, exponential scheduling
	alpha_upd = max(0.5, math.exp(-0.7*T/max(sigma, 1.e-4)))
	Tnew = alpha_upd*T
	return Tnew
  
SEED = 2
np.random.seed(SEED)

TotalCalls = 0

curTrials = 0.0
curAccept = 0.0

T = getInitTemp(); print 'Initial temp:', T
x = x_init()
D = D_init()

T_Values = [objFunction(x)] #Storing values accepted for current temperature
AllValues = [] 				#Storing all the values we accept, for later analysis
AllSteps = []				#Storing the magnitudes of the performed steps
PlotSteps = []; PlotSteps.append(x) #Storing visited locations

while (TotalCalls <= CallsLimit):

  if (curTrials > TrialsNumber or curAccept > AcceptsNumber): #The condition to update the temperature
    sigma = np.std(T_Values)
    T = T_upd(T, sigma)
    acceptRate = curAccept/curTrials*100.0 #Acceptence rate with this temperature
    curAccept = 0; curTrials = 0		   
    T_Values = [T_Values[-1]]			   #Storing the last value we got on previous temperature cycle

  curValue = T_Values[-1]

  dx = MakeStep(x, D)
  xNew = [x[i] + dx[i] for i in range(dimension)]
  newValue = objFunction(xNew)

  stepSize = np.linalg.norm(dx)
  if (newValue > curValue):
    AcceptProbability = math.exp(-(newValue-curValue)/(T*stepSize)) #Calculating step size, temperature and df probability of acceptance, if new position is worse
  else: 
    AcceptProbability = 1.0


  if (np.random.rand() < AcceptProbability): #If we accept the step
    T_Values.append(newValue) 			     #Appending obj func value
    x = xNew
    curAccept+=1.0
    D = D_upd(D, dx) 						#Updating D matrix
    AllValues.append(newValue)				#Storing
    AllSteps.append(stepSize)
    PlotSteps.append(x)
  curTrials += 1.0							


print 'Best point: ', PlotSteps[np.argmin(AllValues)+1]
print 'Best value: ', np.min(AllValues)