added schaffer2 function

added output for improvement of fitness
added CMAES algorithm, that requires further testing

functions.py

@@ -17,17 +17,29 @@ class Function:
         X, Y = np.meshgrid(x, y)
         return x, y, self.eval(X, Y)
 
+
 class Rastrigin(Function):
     def __init__(self, A: int = 10):
         super().__init__(
-                xlim=(-5.12, 5.12), 
+                xlim=(-6.12, 6.12), 
-                ylim=(-5.12, 5.12), 
+                ylim=(-6.12, 6.12), 
                 minimum=(0, 0),
                 eval=lambda x, y: self.A * 2 + \
                         (x**2 - self.A * np.cos(2 * np.pi * x)) + \
                         (y**2 - self.A * np.cos(2 * np.pi * y)))
         self.A = A
 
+
+class Schaffer2(Function):
+    def __init__(self):
+        super().__init__(
+                xlim=(-50, 50),
+                ylim=(-50, 50),
+                minimum=(0, 0),
+                eval=lambda x, y: 0.5 + \
+                        (np.sin(x**2 - y**2)**2 - 0.5) / \
+                        ((1 + 0.001 * (x**2 + y**2))**2))
+
 class Sphere(Function):
     def __init__(self):
         super().__init__(

main.py

@@ -0,0 +1,43 @@
+import functions
+import numpy as np
+import solvers
+
+from matplotlib import pyplot as plt
+
+
+def main():
+    plot_rest = True
+    f = functions.Rastrigin()
+    # f = functions.Sphere()
+    # s = solvers.SimpleEvolutionStrategy(f, np.array([2, 2]))
+    s = solvers.CMAEvolutionStrategy(f, np.array([2.0, 2.0]), np.array([[2.0, 0.0], [0.0, 2.0]]))
+    old_fitness = 100
+    fitness = old_fitness * 0.9
+    old = None
+    plt.set_cmap('Spectral')
+
+    while abs(old_fitness - fitness) > 0.001:
+        old_fitness = fitness
+        samples = s.sample(2000)
+        elite, fitness = s.rank(samples, 150)
+        s.update(elite)
+
+        if plot_rest:
+            rest = np.setdiff1d(samples, elite, assume_unique=True)
+            rest = rest.reshape((int(rest.shape[0]/2), 2))
+        plt.pcolormesh(*f.grid())
+        if plot_rest:
+            plt.scatter(*rest.transpose(), color="dimgrey")
+        if old is not None:
+            plt.scatter(*old.transpose(), color="orange")
+        plt.scatter(*elite.transpose(), color="yellow")
+        # plt.scatter(*elite[0].transpose(), color="green")
+        print('old fitness: {}\nnew fitness: {}\nimprovement: {}\n'.format(old_fitness, fitness,
+            old_fitness - fitness))
+        plt.show()
+        old = elite 
+
+
+if __name__ == '__main__':
+    main()
+

solvers.py

@@ -10,7 +10,7 @@ class Solver(ABC):
         self.function = function
 
     @abstractmethod
-    def rank(self, samples: Iterable[Iterable[float]], elite_size: int) -> (Iterable[Tuple], Iterable[Tuple], float):
+    def rank(self, samples: Iterable[Iterable[float]], elite_size: int) -> (Iterable[Iterable[float]], float):
         pass
 
     @abstractmethod
@@ -18,7 +18,7 @@ class Solver(ABC):
         pass
 
     @abstractmethod
-    def update(elite: Iterable[Tuple]):
+    def update(elite: Iterable[Iterable[float]]):
         pass
 
 
@@ -33,7 +33,7 @@ class SimpleEvolutionStrategy(Solver):
         else:
             self.sigma = [1] * len(mu)
 
-    def rank(self, samples: Iterable[Iterable[float]], elite_size: int) -> (Iterable[Tuple], Iterable[float]):
+    def rank(self, samples: Iterable[Iterable[float]], elite_size: int) -> (Iterable[Iterable[float]], Iterable[float]):
         fitness = self.function.eval(*samples.transpose())
         samples = samples[np.argsort(fitness)]
         fitness = np.sort(fitness)
@@ -43,6 +43,41 @@ class SimpleEvolutionStrategy(Solver):
     def sample(self, n: int) -> Iterable[Iterable[float]]:
         return np.array([np.random.multivariate_normal(self.mu, np.diag(self.sigma)) for _ in range(n)])
 
-    def update(self, elite: Iterable[Tuple]):
+    def update(self, elite: Iterable[Iterable[float]]):
         self.mu = elite[0]
 
+
+class CMAEvolutionStrategy(Solver):
+    def __init__(self, function: Function, mu: Iterable[float], covariances: Iterable[float] = None):
+        if covariances is not None and len(mu) != covariances.shape[0] and len(mu) != covariances.shape[1]:
+            raise Exception('Length of mu and covariance matrix must match')
+        super().__init__(function)
+        self.mu = mu
+        if covariances is not None:
+            self.covariances = covariances
+        else:
+            self.covariances = np.array([[1.0, 0.0], [0.0, 1.0]])
+
+    def rank(self, samples: Iterable[Iterable[float]], elite_size: int) -> (Iterable[Iterable[float]], Iterable[float]):
+        fitness = self.function.eval(*samples.transpose())
+        samples = samples[np.argsort(fitness)]
+        fitness = np.sort(fitness)
+        elite = samples[0:elite_size]
+        return elite, fitness[0]
+
+    def sample(self, n: int) -> Iterable[Iterable[float]]:
+        return np.array([np.random.multivariate_normal(self.mu, self.covariances) for _ in range(n)])
+
+    def update(self, elite: Iterable[Iterable[float]]):
+        mu = self.mu
+        x = elite.transpose()[0]
+        y = elite.transpose()[1]
+        self.mu[0] = np.average(x)
+        self.mu[1] = np.average(y)
+        
+        # TODO fix covariance matrix calculation using the tutorial by Nikolaus Hansen
+        self.covariances[0][0] = np.average((x - mu[0])**2)
+        self.covariances[1][1] = np.average((y - mu[1])**2)
+        self.covariances[0][1] = np.average((x - mu[0]) * (y - mu[1]))
+        # self.covariances = np.cov(elite.transpose())
+

test.py

@@ -1,35 +0,0 @@
-import functions
-import numpy as np
-import solvers
-
-from matplotlib import pyplot as plt
-
-
-plot_rest = True
-f = functions.Rastrigin()
-# f = functions.Sphere()
-s = solvers.SimpleEvolutionStrategy(f, np.array([2, 2]))
-old_fitness = 100
-fitness = old_fitness * 0.9
-old = None
-plt.plasma()
-
-while abs(old_fitness - fitness) > 0.001:
-    old_fitness = fitness
-    samples = s.sample(100)
-    elite, fitness = s.rank(samples, 10)
-    s.update(elite)
-
-    if plot_rest:
-        rest = np.setdiff1d(samples, elite, assume_unique=True)
-        rest = rest.reshape((int(rest.shape[0]/2), 2))
-    plt.pcolormesh(*f.grid())
-    if plot_rest:
-        plt.scatter(*rest.transpose(), color="dimgrey")
-    if old is not None:
-        plt.scatter(*old.transpose(), color="lightgrey")
-    plt.scatter(*elite.transpose(), color="yellow")
-    plt.scatter(*elite[0].transpose(), color="green")
-    plt.show()
-    old = elite 
-