Genetic Algorithm over Complex Geometries — Torus, Klein-Bottle, Bumpy Sphere and more!
Uncover the beautiful dance between optimisation theory and topology
This tutorial will teach you how to build your own Genetic Algorithm (GA) from scratch using python. To make the things more interesting, we will solve the problem over a torus (donut shaped manifold). We will also sketch how to apply the GA to other geometries: Klein-Bottles, Sea Shells, Bumpy-spheres, etc.
For a an introduction to the Genetic Algorithm please visit: The Genetic Algorithm Made Simple
Setting the problem
The Genetic Algorithm (GA) works by creating a population of agents that explores a space of parameters until they converge towards a point with maximum fitness (or minimum loss). We call this point the target. We will make this point unique. In contrast to deep reinforcement learning we won’t require the loss/fitness function to be differentiable.
We will create a population of GA agents living on our torus. They will start with random positions over a coordinate patch that we define. Their coordinates (θ, φ) in the the torus will be our update parameters. This parameters will be updated in each iteration creating a sequence {(θ1, φ1), (θ2, φ2), (θ3, φ3), …} that is optimised to approach the target (θ*, φ*) at the torus.
In this way, after each iteration, our population of agents will move closer and closer to the target point.
Topologically, the torus is the cartesian product of two circles: S1 x S1. Thus only two circular coordinates (θ, φ) are necessary to fully characterise any point in this surface. Note all toruses have two radiuses.
Getting Started
Lets start our project with some basic imports
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy as sc
from more_itertools import sort_together
from numpy import sin, cos, exp, pinow lets define the key parameters
num_population = 80
niter = 70
TwoPi = 2*pi
lim_init = [[0, 0.25*TwoPi], [0, 0.25*TwoPi]]
target = [3*pi/2, 2*pi]
objective_function = lambda x: (x[0]%TwoPi- target[0])**2 + (x[1]%TwoPi - target[1])**2Here we have defined an initial population of 80 agents and a total of 70 GA iterations. The array lim_init contains the starting coordinate patch, target the goal (θ*, φ*) coordinates and objective_function the loss we want GA to minimise. We choose an L2 loss. Note the coordinates x[0], x[1] are made cyclical. This is important. These are the two circles that define a torus. This choice bounds the entire problem to a torus!
Now lets define a parent class to describe our GA agents
class Agent:
""" Agents for the genetic algorithm """
def __init__(self, lim_init=lim_init):
self.x = np.array(
[np.random.uniform(lim_init[0][0], lim_init[0][1], 1)[0],
np.random.uniform(lim_init[1][0], lim_init[1][1], 1)[0]])
self.selected = False
self.fitness = 0
def calculate_fitness(self, objective_function):
self.fitness = 1.0 - objective_function(self.x)
def clear(self):
self.selected = False
def print_state(self):
print(f'x = {self.x}, selected={self.selected}')our agents will be instances of this parent class, thus will have methods for initialisation (in a random state), calculate_fitness and print_state. The class properties are just:
- x: to hold the actual coordinates / optimisation parameters
- selected: boolean indicating if the current agent is selected for mating
- fitness: to hold the fitness value associate with the coordinates x.
The GA goal is to maximise the fitness.
Building the Genetic Algorithm
The schematic procedure behind the GA can be found in the figure below
For the GA to work we need to construct 3 important functions: select, mate and mutate. Then we will build an iterative pipeline using this building blocks.
Lets dive into the necessary functions
Select
def select(population: list, frac_top=0.8, frac_bottom=0.1) -> list:
""" frac_top: fraction of the population selected for breeding
frac_bottom: : Breeding bad ones for diversity
"""
assert 0.0 <= frac_top <= 1.0
assert 0.0 <= frac_top <= 1.0
assert 0.0 <= frac_top+frac_bottom <= 1.0
# calculate fitness
for i, _ in enumerate(population):
population[i].calculate_fitness(objective_function)
# sort according to increasing fitness
population_fitness = [ pop.fitness for pop in population ]
population = sort_together([population_fitness, population])[1]
# select the top frac elements for breeding
cutoff_nr_top = int(frac_top*len(population))
for i, _ in enumerate(population):
if i<= cutoff_nr_top: continue
population[i].selected = True
# select also a random sample of the other elements for breeding
cutoff_nr_bottom = int(frac_bottom*len(population))
for i in range(cutoff_nr_bottom):
random_index = np.random.randint(0, cutoff_nr_top-1)
population[random_index].selected = True
return populationStart by calculation the fitness over the entire population. Then sort the agents based on their fitness. Select the top X% of the population for breeding. For parameter diversity we keep also a few samples from the the bottom population with lower fitness scores. The property selected = True labels the agents that will go into mating.
Mate
def cross_over(pop_old: list) -> list:
pop_new = []
population_nr = 0
while population_nr < len(pop_old):
# Random indices
i_random = np.random.randint(0, len(pop_old))
j_random = np.random.randint(0, len(pop_old))
# Skip agents that were not selected
if (pop_old[i_random].selected is False) or (pop_old[j_random].selected is False): continue
# Breed random pairs of agents
child = breed(pop_old[i_random], pop_old[j_random])
# Add new agents to the pop_new
pop_new.append(child)
population_nr += 1
return pop_new
def breed(parent1: Agent, parent2: Agent) -> Agent:
child = Agent()
child.x[0] = np.random.choice([parent1.x[0], parent2.x[0]])
child.x[1] = np.random.choice([parent1.x[1], parent2.x[1]])
return childNow we take random pairs of selected agents and we breed them together. The breeding function creates new coordinates x[0], x[1] by picking the values from the parents. This choice is made randomly. Ex: breed([1,1], [2,2]) = [1,2]. Thus this function is non-deterministic. We breed random pairs until we have matched the desired number of new agents.
Mutate
def mutate(ind: Agent, mutation_proba=0.75, mutation_strenght=0.31) -> Agent:
if np.random.rand()>=mutation_proba:
for i, _ in enumerate(ind.x):
ind.x[i] = ind.x[i] + np.random.uniform(-mutation_strenght, mutation_strenght)
return indMutation is simple. Just add a random number mutation_strenght (positive or negative) to the current coordinates. In our case we will pick this random number from a uniform distribution. Not all agents need to be mutated, mutation_proba will decide randomly who gets to be mutated.
The stronger the mutations, the more “explorative” the agents become. It means that the child population will deviate more from the parent population. This is desired in cases where several fitness maxima are present, and we wish to converge to the biggest fitness maxima. These type of issues are often the pain-point of “local” optimisation methods like gradient descent.
Main Pipeline — Genetic Algorithm
# initialize population
pop_new = [Agent() for i in range(num_population)]
# initialize the snapshots (for the movie)
snapshot = []
# Run iterations
for t in range(niter):
# replacing old population with new one
pop_old = pop_new
# selection
pop_old = select(pop_old)
# breeding (done by selecting pais at random)
pop_new = cross_over(pop_old)
# mutation
for i, _ in enumerate(pop_new):
pop_new[i] = mutate(pop_new[i])
# clear state
for i, _ in enumerate(pop_new):
pop_new[i].clear()
# save snapshot coordinates
snapshot.append([pop.x for pop in pop_new])Time to run our GA pipeline. We will save snapshots of our population in the snapshot array for plotting purposes later. After the pipeline runs, our GA finds a population of agents with coordinates {(θ, φ)} that should be close to the target point (θ*, φ*) in the torus.
Extracting the best result
# Calculate fitness
for i, _ in enumerate(pop_new):
pop_new[i].calculate_fitness(objective_function)
population_fitness = [ pop.fitness for pop in pop_new ]
# Extract best fitness
index_best = np.argmax(population_fitness)
# print the best one!
print('best solution coordinates:')
pop_new[index_best].print_state()
print('best solution fitness:')
print('x_fitness = ', pop_new[index_best].fitness)The best agents (the ones closes to the target) can be now found in our end population. Is just a matter of calculating the fitness of all agents and picking the one with higher fitness!
Find below the results I obtained. As the starting points are random, each run is unique.
Let’s jump to 3D!
Now that we have obtained the results of our model, is nice to see visualise our results in 3D. To do this we bring the intrinsic manifold coordinates of the torus (u, v) = (θ, φ) to 3D cartesian (x, y, z) coordinates.
The function torus, will allow such conversion:
def torus(u, v, r=0.5, R=1.):
X = (R + r * cos(v)) * cos(u)
Y = (R + r * cos(v)) * sin(u)
Z = r * np.sin(v)
return X, Y, ZNote that this is the first time that we have to include the two radiuses of the torus: (r, R). We had no need of these numbers in the GA calculation before. This is because all toruses are topologically equivalent, independently of their radiuses.
Population migration towards the target
It’s time to bring all together with a nice animation showing the GA population migrating from their starting positions towards the target point.
Let’s use our “snapshots” array to plot the agent coordinates over the toroidal surface. Define
def set_axis(ax):
x_min, x_max, y_min, y_max, z_min, z_max = (-1,1,-1,1,-1,1)
ax.set_xlim3d(x_min, x_max); ax.set_ylim3d(y_min, y_max); ax.set_zlim3d(z_min, z_max)to get nicely scaled plots. Then we proceed on plotting with
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib.animation
import numpy as np
from IPython.display import display, clear_output
# define figure
fig, ax = plt.subplots(figsize=(12,12))
ax = fig.gca(projection = '3d')
set_axis(ax)
# define parametric surface
parameter = np.linspace(0, 2*TwoPi, 132)
u, v = np.meshgrid(parameter, parameter)
X, Y, Z = torus(u, v)
# Plot animation based on snapshots
for it in range(niter):
ax.clear(); set_axis(ax)
# generate population & target data
u = np.array([ pop_x[0] for pop_x in snapshot[it]])
v = np.array([ pop_x[1] for pop_x in snapshot[it]])
# map population intrinsecal coordinates to 3D embedding coordinates
x, y, z = torus(u, v)
tx, ty, tz = torus(target[0], target[1])
# plot 3D
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, color='w', linewidth=0, antialiased=False, alpha=0.15)
ax.scatter(xs=x, ys=y, zs=z, s=70, c='k')
ax.scatter(xs=tx, ys=ty, zs=tz, s=200, c='r')
# clear when the next frame is ready
clear_output(wait=True)
# display current
display(fig)the result should look like this:
If we plot different generations of agents with different colors we will get the plot below. This shows the population migration towards the target:
More manifolds…
If you followed this tutorial carefully, you can now try implementing some other examples yourself, here I show you some examples using the Klein-Bottle, Bumpy-Sphere, Sea-Shell, and the Paraboloid
You can reproduce the 3D embedding of these manifolds with the following python functions:
def paraboloid(X, Y):
Z = X**2 + Y**2
return X, Y, Z
def seashell(u, v):
X = 2*(1-exp(u/(6*pi)))*cos(u)*cos(v/2)**2;
Y = 2*(-1+exp(u/(6*pi)))*sin(u)*cos(v/2)**2;
Z = 1-exp(u/(3*pi))-sin(v)+exp(u/(6*pi))*sin(v);
return X, Y, Z
def bumpy_sphere(u, v):
R = 1 + 0.1*sin(9*u)*cos(9*v)
X = R*cos(u)*sin(v)
Y = R*sin(u)*sin(v)
Z = R*cos(v)
return X, Y, Z
def klein_bottle(u, v):
X = -(2/15)*cos(u)*(3*cos(v)-30*sin(u)+90*cos(u)**4*sin(u)- \
60*cos(u)**6*sin(u)+5*cos(u)*cos(v)*sin(u));
Y = -(1/15)*sin(u)*(3*cos(v)-3*cos(u)**2*cos(v)-48*cos(u)**4*cos(v)+ \
48*cos(u)**6*cos(v)-60*sin(u)+5*cos(u)*cos(v)*sin(u)
-5*cos(u)**3*cos(v)*sin(u) -80*cos(u)**5*cos(v)*sin(u)+ \
80*cos(u)**7*cos(v)*sin(u));
Z = (2/15)*(3+5*cos(u)*sin(u))*sin(v);
return X, Y, ZPlease try this yourself. Take into account that you need to modify the objective function to match the topology of the manifold you want to reproduce. Keep in mind which coordinates are linear and which are cyclical (Hint: The fundamental polygon concept might help: https://en.wikipedia.org/wiki/Fundamental_polygon).
Good luck and keep on learning!
