initiation list04

list04/functions.py

@@ -0,0 +1,37 @@
+import numpy as np
+
+
+def objective_function_F1(X):
+    # Sphere function (minimum at 0)
+    return - np.sum(X**2, axis=1)
+
+
+def objective_function_F1a(X):
+    # Sphere function - modified
+    return - (X[:, 0]**2 + 9*X[:, 1]**2)
+
+
+def objective_function_F1b(X):
+    # Sphere function - modified
+    return - (X[:, 0]**2 + 625*X[:, 1]**2)
+
+
+def objective_function_F1c(X):
+    # Sphere function - modified
+    return - (X[:, 0]**2 + 2*X[:, 1]**2 - 2 * X[:, 0] * X[:, 1])
+
+
+def objective_function_F6(X):
+    # Rastrigin function (minimum at 0)
+    return - 10.0 * X.shape[1] - np.sum(X**2, axis=1) + 10.0 * np.sum(np.cos(2 * np.pi * X), axis=1)
+
+
+def objective_function_F7(X):
+    # Schwefel function (minimum at 420.9687)
+    # (REMARK: should be considered only on [-500, 500]^d, because there are better minima outside)
+    return - 418.9829 * X.shape[1] + np.sum(X * np.sin(np.sqrt(np.abs(X))), axis=1)
+
+
+def objective_function_F8(X):
+    # Griewank function (minimum at 0)
+    return - 1 - np.sum(X**2 / 4000, axis=1) + np.prod(np.cos(X / np.sqrt(np.linspace(1, X.shape[1], X.shape[1]))), axis=1)

list04/mutations.py

@@ -0,0 +1,75 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
+from functions import *
+
+
+def plot_mutation(
+        mutations,
+        domain_X=np.arange(-5, 5, 0.25),
+        domain_Y=np.arange(-5, 5, 0.25),
+        objective_function=objective_function_F1a,
+        original_individual=np.array([[1, 1]]),
+):
+    X, Y = np.meshgrid(domain_X, domain_Y)
+    Z = - objective_function(np.vstack([X.ravel(),
+                                        Y.ravel()]).T).reshape(X.shape[0], X.shape[1])
+    plt.figure(figsize=(6, 6))
+    plt.contour(X, Y, Z, 50)
+    plt.plot(mutations[:, 0], mutations[:, 1], 'ro')
+    plt.plot(original_individual[0, 0],
+             original_individual[0, 1], 'k*', markersize=24)
+    # plt.title('Mutation 1')
+    plt.show()
+
+
+def mutate(
+    sigma,
+    original_individual=np.array([[1, 1]]),
+    N=250,
+    d=2,
+):
+    if (
+            isinstance(sigma, int)
+            or isinstance(sigma, float)
+            or len(sigma.shape) == 1
+    ):
+        return original_individual + sigma * np.random.randn(N, d)
+
+    return (
+        original_individual
+        + np.dot(
+            np.random.randn(N, d),
+            np.linalg.cholesky(sigma).T
+        )
+    )
+
+
+def benchmark_mutation(
+    sigma,
+    objective_function=objective_function_F1a,
+    original_individual=np.array([[1, 1]])
+):
+    indv_res = objective_function(original_individual)[0]
+
+    x = np.array([objective_function(mutate(sigma)) for i in range(100)])
+    better_than_indv = (x < indv_res).sum(axis=1)
+    best_from_pop = x.max(axis=1)
+
+    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3))
+
+    ax1.hist(better_than_indv / 250)
+    ax1.set_title('Amount of better children')
+    ax2.hist(best_from_pop)
+    ax2.set_title('Bests From Children')
+    fig.suptitle(
+        objective_function.__name__,
+        verticalalignment='top',
+        fontsize=15,
+        y=1.05
+
+    )
+    fig.tight_layout()
+    plt.show()

list04/notebook_examples.py

@@ -0,0 +1,34 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
+from notebook_utils import *
+from functions import *
+
+N = 250
+d = 2
+
+objective_function = objective_function_F1a
+original_individual = np.array([[1, 1]])
+
+
+def example_1():
+    sigma = 0.25
+    mutations = original_individual + sigma * np.random.randn(N, d)
+    domain_X = np.arange(-5, 5, 0.25)
+    domain_Y = np.arange(-5, 5, 0.25)
+    X, Y = np.meshgrid(domain_X, domain_Y)
+    Z = - objective_function(
+        np.vstack([
+            X.ravel(),
+            Y.ravel()
+        ]).T
+    ).reshape(X.shape[0], X.shape[1])
+
+    plt.figure(figsize=(9, 9))
+    plt.contour(X, Y, Z, 50)
+    plt.plot(mutations[:, 0], mutations[:, 1], 'ro')
+    plt.plot(original_individual[0, 0],
+             original_individual[0, 1], 'k*', markersize=24)
+    plt.title('Mutation 1')
+    plt.show()

list04/notebook_utils.py

@@ -0,0 +1,211 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
+from functions import *
+
+
+def es(
+    objective_function,
+    chromosome_length,
+    population_size,
+    number_of_iterations,
+    number_of_offspring,
+    number_of_parents,
+    sigma,
+    tau,
+    tau_0,
+    log_frequency=1,
+    verbose=True
+):
+
+    best_solution = np.empty((1, chromosome_length))
+    best_solution_objective_value = 0.00
+
+    log_objective_values = np.empty((number_of_iterations, 4))
+    log_best_solutions = np.empty((number_of_iterations, chromosome_length))
+    log_best_sigmas = np.empty((number_of_iterations, chromosome_length))
+
+    # generating an initial population
+    current_population_solutions = 100.0 * \
+        np.random.rand(population_size, chromosome_length)
+    current_population_sigmas = sigma * \
+        np.ones((population_size, chromosome_length))
+
+    # evaluating the objective function on the current population
+    current_population_objective_values = objective_function(
+        current_population_solutions)
+
+    for t in range(number_of_iterations):
+
+        # selecting the parent indices by the roulette wheel method
+        fitness_values = current_population_objective_values - \
+            current_population_objective_values.min()
+        if fitness_values.sum() > 0:
+            fitness_values = fitness_values / fitness_values.sum()
+        else:
+            fitness_values = 1.0 / population_size * np.ones(population_size)
+        parent_indices = np.random.choice(
+            population_size,
+            (number_of_offspring, number_of_parents),
+            True,
+            fitness_values
+        ).astype(np.int64)
+
+        # creating the children population by Global Intermediere Recombination
+        children_population_solutions = np.zeros(
+            (number_of_offspring, chromosome_length)
+        )
+        children_population_sigmas = np.zeros(
+            (number_of_offspring, chromosome_length)
+        )
+        for i in range(number_of_offspring):
+            children_population_solutions[i, :] = current_population_solutions[parent_indices[i, :], :].mean(
+                axis=0)
+            children_population_sigmas[i, :] = current_population_sigmas[parent_indices[i, :], :].mean(
+                axis=0)
+
+        # mutating the children population by adding random gaussian noise
+        children_population_sigmas = (
+            children_population_sigmas
+            * np.exp(
+                tau * np.random.randn(number_of_offspring, chromosome_length)
+                + tau_0 * np.random.randn(number_of_offspring, 1)
+            )
+        )
+        children_population_solutions = children_population_solutions + \
+            children_population_sigmas * \
+            np.random.randn(number_of_offspring, chromosome_length)
+
+        # evaluating the objective function on the children population
+        children_population_objective_values = objective_function(
+            children_population_solutions
+        )
+
+        # replacing the current population by (Mu + Lambda) Replacement
+        current_population_objective_values = np.hstack(
+            [
+                current_population_objective_values,
+                children_population_objective_values
+            ]
+        )
+        current_population_solutions = np.vstack(
+            [current_population_solutions, children_population_solutions])
+        current_population_sigmas = np.vstack(
+            [current_population_sigmas, children_population_sigmas])
+
+        I = np.argsort(current_population_objective_values)[::-1]
+        current_population_solutions = current_population_solutions[I[:population_size], :]
+        current_population_sigmas = current_population_sigmas[I[:population_size], :]
+        current_population_objective_values = current_population_objective_values[
+            I[:population_size]]
+
+        # recording some statistics
+        if best_solution_objective_value < current_population_objective_values[0]:
+            best_solution = current_population_solutions[0, :]
+            best_solution_objective_value = current_population_objective_values[0]
+        log_objective_values[t, :] = [
+            current_population_objective_values.min(),
+            current_population_objective_values.max(),
+            current_population_objective_values.mean(),
+            current_population_objective_values.std()
+        ]
+        log_best_solutions[t, :] = current_population_solutions[0, :]
+        log_best_sigmas[t, :] = current_population_sigmas[0, :]
+
+        if verbose and np.mod(t + 1, log_frequency) == 0:
+            print("Iteration %04d : best score = %0.8f, mean score = %0.8f." % (
+                t + 1, log_objective_values[:t+1, 1].max(), log_objective_values[t, 2]))
+
+    return (
+        best_solution_objective_value,
+        best_solution,
+        log_objective_values,
+        log_best_solutions,
+        log_best_sigmas
+    )
+
+
+def plot_es(res):
+    (
+        best_objective_value,
+        best_chromosome,
+        history_objective_values,
+        history_best_chromosome,
+        history_best_sigmas
+    ) = res
+    plt.figure(figsize=(18, 4))
+    plt.plot(history_objective_values[:, 0], label='min')
+    plt.plot(history_objective_values[:, 1], label='max')
+    plt.plot(history_objective_values[:, 2], label='mean')
+    plt.xlabel('iteration')
+    plt.ylabel('objective function value')
+    plt.title('min/avg/max objective function values')
+    plt.legend()
+    plt.show()
+
+    plt.figure(figsize=(18, 4))
+    plt.plot(history_best_sigmas)
+    plt.xlabel('iteration')
+    plt.ylabel('sigma value')
+    plt.title('best sigmas')
+    plt.show()
+
+
+def run_es(
+    d,
+    N,
+    T,
+    func,
+    log_frequency=10,
+    number_of_parents=2,
+    plot=True,
+    sigma=50.0,
+    verbose=True,
+
+):
+    result = es(
+        objective_function=func,
+        chromosome_length=d,
+        population_size=N,
+        number_of_iterations=T,
+        number_of_offspring=2*N,
+        number_of_parents=number_of_parents,
+        sigma=sigma,
+        tau=1/np.sqrt(2*d),
+        tau_0=1/np.sqrt(2*np.sqrt(d)),
+        log_frequency=log_frequency,
+        verbose=verbose
+    )
+    (
+        best_objective_value,
+        best_chromosome,
+        history_objective_values,
+        history_best_chromosome,
+        history_best_sigmas
+    ) = result
+    if plot:
+        plot_es(result)
+    return result
+
+
+def plot_3D_benchmark_function(objective_function, domain_X, domain_Y, title):
+    plt.figure(figsize=(12, 8))
+    ax = plt.gca(projection='3d')
+    X, Y = np.meshgrid(domain_X, domain_Y)
+    Z = - objective_function(np.vstack([X.ravel(),
+                                        Y.ravel()]).T).reshape(X.shape[0], X.shape[1])
+    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1,
+                           cmap=cm.hot, linewidth=0, antialiased=True)
+    plt.title(title)
+    plt.show()
+
+
+def plot_contour_benchmark_function(objective_function, domain_X, domain_Y, title):
+    plt.figure(figsize=(9, 9))
+    X, Y = np.meshgrid(domain_X, domain_Y)
+    Z = - objective_function(np.vstack([X.ravel(),
+                                        Y.ravel()]).T).reshape(X.shape[0], X.shape[1])
+    plt.contour(X, Y, Z, 50)
+    plt.title(title)
+    plt.show()