#!/usr/bin/env python

"""
This is a simple Python implementation of the (mu/mu_I, lambda)-sigmaSA-ES
as discussed in 
http://www.scholarpedia.org/article/Evolution_Strategies
Based on
https://homepages.fhv.at/hgb/downloads/mu_mu_I_lambda-ES.oct
"""

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


MU = 3                    # number of parents
LAMBDA = 12               # number of offspring

X_INIT = np.ones(1)       # initial parent vector 
#X_INIT = np.ones(2)       # initial parent vector 
SIGMA_INIT = 1            # initial global mutation strength sigma 

SIGMA_MIN = 1e-10         # ES stops when sigma is smaller than sigma_min

###########################################################

# Function to be optimized (sphere test function as an example)
# individual: the individual to evaluate
def sphere_function(indiv):
    x = indiv.x
    assert x.shape == X_INIT.shape, x
    y = np.dot(x, x)
    return y

fitness = sphere_function

###########################################################

# The individual class
# x: the individual's value
# sigma: the individual's sigma
class Individual():
    def __init__(self, x, sigma):
        self.x = x
        self.sigma = sigma
        self.cost = fitness(self)

    def __str__(self):
        return "{0} {1} {2}".format(self.x, self.sigma, self.cost)

###########################################################

# This sorts the population according to the individuals' fitnesses
# pop: a list of Individual objects
def select_individuals(pop):
    pop.sort(key=lambda indiv: indiv.cost, reverse=False)
    return pop[:MU]

# This performs intermediate (multi-) recombination
# parents: a list of Individual objects
def recombine_individuals(parents):
    parents_y = np.array([indiv.x for indiv in parents])
    parents_sigma = np.array([indiv.sigma for indiv in parents])
    recombinant = Individual(parents_y.mean(axis=0), parents_sigma.mean())
    return recombinant

###########################################################

def main():

    # Initialization
    n = X_INIT.shape[0]         # determine search space dimensionality n   
    tau = 1. / math.sqrt(2.*n)  # self-adaptation learning rate

    # Initializing individual population
    parent_pop = [Individual(X_INIT, SIGMA_INIT) for i in range(MU)]

    # Evolution loop of the (mu/mu_I, lambda)-sigma-SA-ES
    while parent_pop[0].sigma > SIGMA_MIN:
        offspring_pop = []
        recombinant = recombine_individuals(parent_pop) # TODO: BUG ? this statement may be in the next line
        for i in range(1, LAMBDA):
            offspring_sigma = recombinant.sigma * math.exp(tau * random.normalvariate(0,1))
            offspring_y = recombinant.x + offspring_sigma * np.random.normal(size=n)
            offspring = Individual(offspring_y, offspring_sigma)
            offspring_pop.append(offspring)
        parent_pop = select_individuals(offspring_pop)
        print(parent_pop[0])

# Remark: Final approximation of the optimizer is in "parent_pop[0].x"
#         corresponding fitness is in "parent_pop[0].cost" and the final 
#         mutation strength is in "parent_pop[0].sigma"

if __name__ == "__main__":
    main()

Add a SAES optimizer script.