AI with Python – Genetic Algorithms

This chapter delves into Genetic Algorithms (GAs), a fascinating class of optimization algorithms inspired by Charles Darwin’s theory of natural evolution. We’ll explore their principles and implement them in Python using the DEAP library.

What are Genetic Algorithms?

Genetic Algorithms are search-based heuristics that mimic the process of natural selection and genetics. They belong to the larger family of Evolutionary Algorithms (EAs). GAs are particularly well-suited for solving complex optimization and search problems where traditional methods might struggle.

Developed primarily by John Holland and his collaborators, GAs operate on a population of potential solutions (called individuals or chromosomes). Through mechanisms like selection, crossover (recombination), and mutation, the population evolves over generations, progressively producing better solutions.

Key Concepts in GAs:

Population: A set of candidate solutions to the problem.
Chromosome (Individual): A single candidate solution, typically represented as a string of genes (e.g., binary string, list of numbers).
Gene: A basic building block of a chromosome, representing a part of the solution.
Fitness Function: A function that evaluates how good a solution (chromosome) is. Higher fitness indicates a better solution.
Selection: The process of choosing individuals from the current population to be parents for the next generation. Fitter individuals have a higher chance of being selected (Survival of the Fittest).
Crossover (Recombination): Combines genetic material from two parent chromosomes to create one or more offspring. This allows exploration of new solution regions.
Mutation: Randomly alters one or more genes in a chromosome. This introduces diversity into the population and helps escape local optima.

GAs balance exploration (searching new areas of the solution space) and exploitation (focusing on promising areas found so far).

Using GAs for Optimization Problems

Optimization involves finding the best solution from a set of feasible solutions, typically by maximizing or minimizing an objective function. GAs provide a robust framework for this.

Textual description of an optimization process diagram: An optimization process can be visualized as a loop. It starts with an ‘Initial Guess/Solution’. This solution is fed into a ‘System Model/Objective Function’ which produces ‘Performance Metrics’. These metrics are then analyzed by an ‘Optimization Algorithm’ (like GA). The algorithm generates a ‘New/Improved Solution’, which then feeds back as the new ‘Initial Guess/Solution’, and the loop continues until a satisfactory solution is found or a stopping criterion is met.

General Steps in a Genetic Algorithm

Initialization: Create an initial population of random candidate solutions.
Evaluation: Calculate the fitness of each individual in the population using the fitness function.
Selection: Select parent individuals based on their fitness.
Crossover: Apply crossover operator to pairs of selected parents to create offspring.
Mutation: Apply mutation operator to the offspring (and sometimes to parents) with a certain probability.
Replacement: Form a new population for the next generation, typically by replacing some or all of the old population with the new offspring (and sometimes elite individuals from the previous generation).
Termination: Repeat steps 2-6 until a termination condition is met (e.g., maximum number of generations, fitness threshold reached, no improvement for several generations).

Installing Necessary Packages: DEAP

For implementing GAs in Python, we’ll use DEAP (Distributed Evolutionary Algorithms in Python), a powerful and flexible library for evolutionary computation. Install it via pip:

pip install deap numpy matplotlib

If you’re using Anaconda, you can also try:

conda install -c conda-forge deap

Implementing Solutions using Genetic Algorithms with DEAP

DEAP provides tools to define individuals, populations, genetic operators, and evolutionary algorithms. Let’s explore two examples.

Example 1: The OneMax Problem (Maximizing Ones in a Bitstring)

The OneMax problem is a classic toy problem for GAs. The goal is to evolve a binary string (a list of 0s and 1s) that contains the maximum possible number of 1s. For a string of length N, the optimal solution is a string of N ones.

Import necessary DEAP components and other libraries:

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

from deap import base, creator, tools, algorithms

# Problem parameters
NUM_BITS = 45  # Length of the bitstring
TARGET_SUM = NUM_BITS # We want all ones

# --- 1. Define Fitness and Individual ---
# Create a Maximizing Fitness class (seeks to maximize a single objective)
# weights=(1.0,) means we want to maximize the first (and only) objective value.
creator.create("FitnessMax", base.Fitness, weights=(1.0,))

# Create an Individual class: a list of bits with the FitnessMax property.
creator.create("Individual", list, fitness=creator.FitnessMax)

# --- 2. Initialize Toolbox ---
# The toolbox will store functions for generating individuals and populations,
# and for genetic operators.
toolbox = base.Toolbox()

# Attribute generator: defines how a single gene (bit) is created (0 or 1).
toolbox.register("attr_bool", random.randint, 0, 1)

# Individual generator: defines how an individual (chromosome) is created.
# tools.initRepeat calls toolbox.attr_bool NUM_BITS times to fill the Individual.
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_bool, NUM_BITS)

# Population generator: defines how a population is created (a list of individuals).
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

# --- 3. Define Evaluation (Fitness) Function ---
# This function calculates how 'fit' an individual is.
# For OneMax, fitness is the number of ones in the bitstring.
def evalOneMax(individual):
    return sum(individual), # Return a tuple, as fitness can be multi-objective

# Register the evaluation function in the toolbox.
toolbox.register("evaluate", evalOneMax)

# --- 4. Register Genetic Operators ---
# Crossover operator: how two individuals combine to create offspring.
# tools.cxTwoPoint performs a two-point crossover.
toolbox.register("mate", tools.cxTwoPoint)

# Mutation operator: how an individual is randomly changed.
# tools.mutFlipBit flips bits with a probability indpb (independent probability).
toolbox.register("mutate", tools.mutFlipBit, indpb=0.05) # 5% chance to flip each bit

# Selection operator: how individuals are chosen to be parents.
# tools.selTournament uses tournament selection (tournsize is the number of individuals in each tournament).
toolbox.register("select", tools.selTournament, tournsize=3)

# --- 5. Run the Genetic Algorithm ---
def run_onemax_ga():
    random.seed(42) # For reproducibility

    pop_size = 300
    population = toolbox.population(n=pop_size)

    # Probability of crossover and mutation
    prob_crossing, prob_mutating = 0.7, 0.2 # CXPB, MUTPB
    num_generations = 40

    print('Evolution process starts')

    # Statistics to keep track of during evolution
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", np.mean)
    stats.register("std", np.std)
    stats.register("min", np.min)
    stats.register("max", np.max)

    # Run the GA using a DEAP algorithm (eaSimple is a standard one)
    # halloffame stores the best individuals found during evolution.
    hof = tools.HallOfFame(1) # Store the single best individual
    population, logbook = algorithms.eaSimple(population, toolbox,
                                              cxpb=prob_crossing, mutpb=prob_mutating,
                                              ngen=num_generations, stats=stats,
                                              halloffame=hof, verbose=True)

    print('\n- Evolution ends -')

    best_individual = hof[0]
    print(f'\nBest individual found: {best_individual}')
    print(f'Number of ones: {sum(best_individual)}')
    print(f'Fitness: {best_individual.fitness.values[0]}')

    # Plot statistics
    gen = logbook.select("gen")
    max_fitness_values = logbook.chapters["fitness"].select("max")
    avg_fitness_values = logbook.chapters["fitness"].select("avg")

    plt.figure(figsize=(10,6))
    plt.plot(gen, max_fitness_values, label="Maximum Fitness")
    plt.plot(gen, avg_fitness_values, label="Average Fitness")
    plt.xlabel("Generation")
    plt.ylabel("Fitness")
    plt.title("Fitness over Generations (OneMax Problem)")
    plt.legend()
    plt.grid(True)
    plt.show()

    return population, logbook, hof

if __name__ == "__main__":
    run_onemax_ga()

Expected Output (abbreviated, actual output will show stats for each generation):

Evolution process starts
gen  nevals  avg      std      min  max
0    300     22.41    3.16842  12   31
1    216     25.3033  2.93396  17   33
... (many generations) ...
39   207     44.9     0.3      44   45
40   221     44.9367  0.24362  44   45

- Evolution ends -

Best individual found: [1, 1, ..., 1] (a list of 45 ones)
Number of ones: 45
Fitness: 45.0

The plot will show two lines: one for maximum fitness and one for average fitness per generation. Both should generally increase, with the maximum fitness reaching the target (45 in this case). This demonstrates the GA’s ability to find the optimal solution.

Practical Applications of GAs: Optimization in scheduling, logistics, engineering design, financial modeling, and training machine learning models (e.g., hyperparameter tuning, feature selection).

Example 2: Symbolic Regression

Symbolic Regression is a type of regression analysis that searches the space of mathematical expressions to find the model that best fits a given dataset. Unlike traditional regression that fits parameters to a pre-specified model structure (e.g., linear regression), symbolic regression evolves the structure of the model itself. Genetic Programming (GP), a specialized form of GA, is often used for this.

Goal: Evolve a mathematical function (e.g., y = f(x)) that best fits a set of (x, y) data points. Let’s say our target function is y = x^3 - 2*x^2 + x - 5, and we’ll generate some noisy data around it.

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

from deap import algorithms, base, creator, tools, gp

# --- 1. Define Primitives (Functions and Terminals) ---
# Primitives are the building blocks for the symbolic expressions (trees).
pset = gp.PrimitiveSet("MAIN", 1) # "MAIN" is the name, 1 is arity (number of inputs, e.g., 'x')
pset.addPrimitive(operator.add, 2)       # Addition (2 arguments)
pset.addPrimitive(operator.sub, 2)       # Subtraction (2 arguments)
pset.addPrimitive(operator.mul, 2)       # Multiplication (2 arguments)
def protected_div(left, right):
    try:
        return left / right
    except ZeroDivisionError:
        return 1 # Return a constant if division by zero
pset.addPrimitive(protected_div, 2)    # Protected division
pset.addPrimitive(operator.neg, 1)       # Negation (1 argument)
# pset.addPrimitive(math.sin, 1)         # Sine (optional)
# pset.addPrimitive(math.cos, 1)         # Cosine (optional)

# Terminals: inputs to the functions (e.g., variables like 'x' or constants)
pset.addEphemeralConstant("rand101", lambda: random.uniform(-1, 1)) # Random constant between -1 and 1
pset.renameArguments(ARG0='x') # Rename the input argument to 'x'

# --- 2. Define Fitness and Individual ---
creator.create("FitnessMin", base.Fitness, weights=(-1.0,)) # We want to MINIMIZE the error
creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)

# --- 3. Initialize Toolbox ---
toolbox = base.Toolbox()
# Expression generator: defines how trees are built.
# gp.genHalfAndHalf combines full trees and grow trees for diversity.
toolbox.register("expr", gp.genHalfAndHalf, pset=pset, min_=1, max_=3) # min/max depth of trees

# Individual and Population generators
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.expr)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

# Function to compile the tree into a callable Python function
toolbox.register("compile", gp.compile, pset=pset)

# --- 4. Define Evaluation (Fitness) Function ---
# Generate some sample data to fit
# Target function: y = x^3 - 2*x^2 + x - 5
NUM_POINTS = 50
X_POINTS = np.linspace(-5, 5, NUM_POINTS)
Y_TARGET = X_POINTS**3 - 2*X_POINTS**2 + X_POINTS - 5 + np.random.randn(NUM_POINTS) * 2 # Add some noise

def evalSymbReg(individual, points_x, points_y):
    # Compile the GP tree into a callable function
    func = toolbox.compile(expr=individual)
    # Calculate the Mean Squared Error (MSE) between predicted and target y values
    try:
        sqerrors = [(func(x) - y)**2 for x, y in zip(points_x, points_y)]
        return math.fsum(sqerrors) / len(points_y), # Return as a tuple
    except (OverflowError, ValueError):
        return float('inf'), # Penalize expressions that cause errors

toolbox.register("evaluate", evalSymbReg, points_x=X_POINTS, points_y=Y_TARGET)

# --- 5. Register Genetic Operators (GP specific) ---
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("mate", gp.cxOnePoint) # One-point crossover for trees
toolbox.register("expr_mut", gp.genFull, min_=0, max_=2) # Mutation can generate small subtrees
toolbox.register("mutate", gp.mutUniform, expr=toolbox.expr_mut, pset=pset)

# Decorators to limit tree height (to prevent bloat)
toolbox.decorate("mate", gp.staticLimit(key=operator.attrgetter("height"), max_value=17))
toolbox.decorate("mutate", gp.staticLimit(key=operator.attrgetter("height"), max_value=17))

# --- 6. Run the Genetic Programming Algorithm ---
def run_symbolic_regression_gp():
    random.seed(60)

    pop_size = 500
    population = toolbox.population(n=pop_size)
    hof = tools.HallOfFame(1)

    stats_fit = tools.Statistics(lambda ind: ind.fitness.values)
    stats_size = tools.Statistics(len) # Tree size (number of nodes)
    mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size)
    mstats.register("avg", np.mean)
    mstats.register("std", np.std)
    mstats.register("min", np.min)
    mstats.register("max", np.max)

    prob_crossing, prob_mutating = 0.8, 0.15
    num_generations = 50

    print('Symbolic Regression GP process starts')
    population, logbook = algorithms.eaSimple(population, toolbox, prob_crossing, prob_mutating,
                                              num_generations, stats=mstats, halloffame=hof,
                                              verbose=True)

    print('\n- GP Evolution ends -')
    best_individual = hof[0]
    print(f'\nBest individual found (tree structure):\n{best_individual}')
    print(f'Fitness (MSE): {best_individual.fitness.values[0]}')

    # Convert the best tree to a string for readability
    best_func_str = str(best_individual)
    print(f'As a string: {best_func_str}')

    # Plot fitness and size over generations
    gen = logbook.select("gen")
    min_fitness_values = logbook.chapters["fitness"].select("min")
    avg_fitness_values = logbook.chapters["fitness"].select("avg")
    avg_size_values = logbook.chapters["size"].select("avg")

    fig, ax1 = plt.subplots(figsize=(12,7))
    line1 = ax1.plot(gen, min_fitness_values, "b-", label="Minimum Fitness (MSE)")
    ax1.plot(gen, avg_fitness_values, "g-", label="Average Fitness (MSE)", alpha=0.5)
    ax1.set_xlabel("Generation")
    ax1.set_ylabel("Fitness (MSE)", color="b")
    ax1.tick_params(axis='y', labelcolor='b')
    ax1.set_yscale('log') # Often useful for fitness

    ax2 = ax1.twinx()
    line2 = ax2.plot(gen, avg_size_values, "r-", label="Average Tree Size")
    ax2.set_ylabel("Average Size", color="r")
    ax2.tick_params(axis='y', labelcolor='r')

    lns = line1 + line2
    labs = [l.get_label() for l in lns]
    ax1.legend(lns, labs, loc="center right")
    plt.title("Symbolic Regression: Fitness and Size over Generations")
    plt.grid(True)
    plt.show()

    # Plot the fit of the best individual
    func_best = toolbox.compile(expr=best_individual)
    Y_PRED = [func_best(x) for x in X_POINTS]
    plt.figure(figsize=(10,6))
    plt.scatter(X_POINTS, Y_TARGET, color='blue', label='Target Data (noisy)')
    plt.plot(X_POINTS, Y_PRED, color='red', linewidth=2, label='GP Evolved Function')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Symbolic Regression: Best Evolved Function vs Target Data')
    plt.legend()
    plt.grid(True)
    plt.show()

    return population, logbook, hof

if __name__ == "__main__":
    run_symbolic_regression_gp()

Expected Output (abbreviated):

Symbolic Regression GP process starts
gen  nevals  avg       std       min       max       avg       std      min  max
0    500     1.11e+06  7.41e+06  12.01     1.01e+08  4.952     1.76492  2    7
1    412     2.31e+05  2.29e+06  9.65      4.64e+07  5.86      2.72186  2    14
... (many generations) ...
49   413     205.133   1071.76   2.0101    2.22e+04  12.604    6.0282   2    31
50   421     24.2023   135.434   2.0095    2901.17   12.55     6.1135   2    30

- GP Evolution ends -

Best individual found (tree structure):
(Example: mul(sub(mul(x, x), x), x))  <-- This is a placeholder, actual tree will be complex
Fitness (MSE): (e.g., 2.0095)
As a string: (e.g., 'mul(sub(mul(x, x), x), x)') <-- Actual string representation of the tree

The first plot will show how fitness (MSE, ideally decreasing) and average tree size change over generations. The second plot will show the original noisy data points and the function evolved by the GP, illustrating how well it fits the data. The evolved function might not be exactly x^3 - 2*x^2 + x - 5 but could be a close approximation or a different functional form that also fits the data well.

Symbolic regression is a powerful tool for automated model discovery in science and engineering.

Further Learning Resources:

DEAP Documentation: https://deap.readthedocs.io/en/master/
A Field Guide to Genetic Programming (by Poli, Langdon, McPhee): A comprehensive resource, often available online.
Introduction to Evolutionary Computing (by Eiben and Smith): A textbook covering GAs and other EAs.