An implementation of Monte Carlo Octree Optimisation

mcoo.py

import math
import numpy as np
import matplotlib.pyplot as plt

# Optimisation functions
# Polynomial
# def f(x):
#   return np.sum(x**2)

# Ackley
def f(x):
  a = 20
  b = 0.2
  c = 2*math.pi
  return -a*np.exp(-b*np.sqrt(np.sum(x**2)/DIMENSIONS)) - np.exp(np.sum(np.cos(c*x))/DIMENSIONS) + a + math.e


# Verbose printing
def print_verbose():
  print x
  print x_fit
  print y
  print y_fit
  print y_hat
  print y_hat_fit
  print min
  print max

# Fitness history graph
def plot_history():
  plt.figure()
  plt.xlabel("Epoch")
  plt.ylabel("Fitness")
  plt.plot(history_fit)
  plt.show()

# Parameters
DIMENSIONS = 2
EPOCHS = 100
DOMAIN_MIN = -5.0
DOMAIN_MAX = 5.0
SEED_FLAG = True
SEED = 1
SWARM_SIZE = 30
VERBOSE = False

if SEED_FLAG:
  np.random.seed(SEED)
  
# Initialize
x = []
x_fit = []
y = None
y_fit = None
y_hat = None
y_hat_fit = None
history_fit = []

min = np.full(DIMENSIONS, DOMAIN_MIN)
max = np.full(DIMENSIONS, DOMAIN_MAX)

for i in range(0, EPOCHS):
  x = []
  x_fit = []

  for j in range(0, SWARM_SIZE):
    # Sample search space
    sample = np.random.random_sample(DIMENSIONS)
    pos = (max - min) * sample + min
    fit = f(pos)
    x.append(pos)
    x_fit.append(fit)

    # Calculate pbest
    if (y is None) or (fit < y_fit):
      y = pos.copy()
      y_fit = fit.copy()

      # Calculate gbest
      if (y_hat is None) or (y_fit < y_hat_fit):
        y_hat = y.copy()
        y_hat_fit = y_fit.copy()

  # Scale search space as octree around gbest
  for k in range(0, DIMENSIONS):
    middle = (min[k] + max[k])/2
    
    # Upper Half
    if (y_hat[k] >= middle):
      min[k] = middle
    else:
      max[k] = middle

  if VERBOSE:
    print_verbose()

  history_fit.append(y_hat_fit)
  print y_hat_fit

plot_history()