I was watching Computerphile's video on using genetic algorithms to solve the Knapsack problem, and I decided to give it a whack.
For anyone running the code, the result_data pandas Data Frame contains the most recent iteration's breakdown of the genomes that are "still alive".
>>> result_data
solution count fitness reproduction_probability population_size
0 11010000 36 11 0.462117 47
1 10000001 1 11 0.462117 1
2 00001100 31 10 0.420788 39
3 01011000 1 10 0.420788 1
4 01000001 1 10 0.420788 1
.. ... ... ... ... ...
90 00100111 1 0 -0.064774 1
91 00011110 1 0 -0.064774 1
92 00010011 1 0 -0.064774 1
93 00001111 1 0 -0.064774 1
94 00001110 1 0 -0.064774 1
The fitness parameter is a particular solution's utility as calculated by the fitness function, which is itself simply the sum of the values of all the items in the solution. When the sum of the weights exceed the specified limit, however, the fitness of that particular solution is set to zero.
I decided the model the probability of successfully reproducing using the hyperbolic tangent function. After each iteration, I calculate the z-score of each solution and I pass that in to the \$\tanh x\$ function, yielding what seemed to me like pretty decent results. Again, I know nothing about this field, though.
I keep track of the extant number of each genome over each iteration, and I use Matplotlib to plot a graph at the end of the genomes which scored higher than an eight on the fitness function (this is pretty arbitrary, but the vast majority of the genomes were massively outdone pretty quickly and became afterthoughts right away, so it helps to not clutter the graph).
I wrote the entire thing as a one-file script, which I know is not ideal, but I was having a blast writing the actual logic. The one change I would make is incorporating configparser and argparse to be able to handle both configuration files and command-line parameters.
The simulation can already handle custom configuration parameters, albeit only like four of them, but the mechanism is really clunky.
I'm excited to hear about any and all potential improvements, even the aforementioned structural refactoring. If I may request one specific area of focus: you might notice that there are what seem like semi-half-hearted attempts at using generators that look like I just gave up and used lists. I haven't really gotten the opportunity to sit down and internalize generators (not to mention coroutines and their async counterparts), so any feedback on that particular subject is particularly appreciated. This project seems like the perfect place for generator expressions to really shine.
Of course, any feedback is always appreciated regardless. Thanks for reading!
"""
Knapsack Problem Genetic Algorithm Solver
"""
from collections import Counter
from dataclasses import dataclass
from itertools import chain, combinations, filterfalse, product, repeat
from operator import attrgetter
from typing import Any, List, Text, TypeVar, Union
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.ticker import MultipleLocator
from more_itertools import unique_everseen
from numpy.random import default_rng
from pandas import DataFrame, Series
def z_scores(data: Series) -> Series:
return (data - data.mean()) / data.std()
@dataclass
class Item:
weight: float = 0.0
value: float = 0.0
class Solution:
def __init__(self, items: List[Item], limit: float, solution: str = None) -> None:
self.items = items
self.limit = limit
self._solution = \
self._generate_random_solution() \
if solution is None else int(solution, base=2)
def __str__(self) -> Text:
return f"{self._solution:08b}"
def __repr__(self) -> Text:
return f"Solution({Text(self)})"
def _generate_random_solution(self) -> int:
global rng
solution = 0
for i, _ in enumerate(self.items):
if rng.binomial(1, 0.5):
solution += (1 << i)
return solution
def fitness(self) -> int:
return solution_fitness(self._solution, self.limit)
def solution_fitness(solution: Text, limit: float) -> float:
# The fitness of a particular population group is equal
# to the total item value that solution would have
# yielded.
#
# Note that that if a particular solution results in
# a weight above the specified limit, the fitness score
# for that solution is automatically graded at zero.
#
fitness = sum(
items[i].value if c == "1" else 0
for i, c in enumerate(f"{solution}"))
# Cap the fitness score of a population group based on
# whether they were able to make it under the allowable
# weight limit. If they weren't, the solution is useless
# and the fitness score for that solution is zero.
return fitness if fitness <= limit else 0.0
@dataclass
class Result:
solution: str
count: int
fitness: float
# This is the maximum weight limit that we are able to carry
# in this scenario. This is the primary constraint against
# all possible solutions will be judged. Any solution which
# goes over this limit (by any amount), will be judged an
# instant failure.
limit = 11.4
# This is the collection from which the solutions will be
# picking. The goal is to come up with the maximum value
# while still coming in under the maximum weight limit,
# is defined above.
items = [
Item(weight=7.0, value=5.0),
Item(weight=2.0, value=4.0),
Item(weight=1.0, value=7.0),
Item(weight=9.0, value=2.0),
Item(weight=6.0, value=4.0),
Item(weight=3.0, value=6.0),
Item(weight=4.0, value=5.0),
Item(weight=5.0, value=6.0)
]
# Define the verbose logging command-line and configuration
# variable.
verbose_logging = False
# Define the number of genetic evolution iterations to
# simulate.
evolution_iterations = 12
# Define the population size.
population_size = 125
# Define the probability of a crossover event occurring.
crossover_probability = 0.57
# Define the probability of random genetic mutations.
genetic_mutation_probability = 0.018
# Define the natural death rate. This is the rate at which
# elements probabilistically die, allowing for weaker
# genetic elements to actually die out.
natural_death_factor = 0.17
# Initialize the random number generator.
rng = default_rng()
# Generate the initial simulated population.
#
# The initial population needs to be randomly generated, but
# every successive iteration thereafter should be created
# based on the genetic evolution of the previous generation.
#
population = (
Solution(items, limit)
for _ in range(population_size))
# Initialize the population variable.
population = list(population)
# Setup pandas series array in which the counts over each
# iteration will be logged for each of the genomes.
series = {
str(solution): Series(str(solution)) for solution in population
}
# Initialize the sentinel value we will use to ensure the
# initial counts get logged only once.
initial_counts_logged = False
# Simulate each successive iteration of population dynamics.
for iteration in range(evolution_iterations):
if verbose_logging:
# Output verbose debugging log indicating that the
# simulation is currently entering iteration n.
print(f"Iteration {iteration:3d}...\n")
# Create a Counter element to determine the relative
# frequency of each possible population element.
counter = Counter(str(element) for element in population)
# Check whether the first counts have been written to
# the pandas series array. If they haven't, go ahead and
# do it.
if initial_counts_logged == False:
# Add the initial iteration's counts to the pandas
# series array.
for key, val in counter.items():
series[key] = series[key].append(Series(val))
# Make sure not to re-log the counts here again,
# since we take care of this at the end of each
# iteration, after performing the genetic mutations
# and cross-over events.
initial_counts_logged = True
# Create a list of the results using the counts mapping
# in the Counter dictionary.
results = [
Result(solution, count, solution_fitness(solution, limit))
for solution, count in counter.items()
]
# Sort the results of the current iteration in reverse
# order, decreasing from highest fitness to lowest
# fitness, then highest count to lowest count, and then
# finally using the solution itself for lexicographic
# sort order.
results = sorted(results, key=attrgetter(
"fitness", "count", "solution"
), reverse=True)
# Create a pandas DataFrame out of the results to make
# statistical analysis way easier.
result_data = DataFrame(results, copy=False)
# Set the probability of successful reproduction for
# each solution by normalizing the fitness of all of the
# solutions in the result data.
#
# At the moment, the standard z-score of the fitness
# for each solution is calculated (the fitness is
# assumed to be Gaussian, although this is absolutely
# not the case), and the result is passed to the
# hyperbolic tangent function.
fitness_z_scores = z_scores(result_data["fitness"])
# In order to transform the z-scores into reproduction
# probabilities, we first add by the absolute value of
# all of the z-scores to ensure we have a range from 0.0
# to max z-score times two. Then, we divide all of the
# z-scores by two times the maximum z-score, in order to
# normalize the values to between 0.0 and 1.0.
#
# Once this has been completed, we finally conclude by
# passing in the normalized values to the hyperbolic
# tangent function, which finally gives us the
# reproduction probabilities for each group.
reproduction_probabilities = \
np.tanh(fitness_z_scores / \
np.max(2.0 * fitness_z_scores))
# Augment the result dataset by adding the reproduction
# probabilities that we just calculated to the data
# frame.
result_data = result_data.assign(
reproduction_probability=reproduction_probabilities)
# I'm just renaming the reproduction probabilities
# series that we calculated above to give it a shorter
# name, in order to comply with Python programming
# standards.
survival_probability = \
result_data["reproduction_probability"]
# Filter out the solutions which have already died out.
#
# This actually causes Python to crash, although I
# haven't taken the time to figure out why, I simply
# disabled this step in the meantime.
#
# result_data = result_data[result_data["reproduction_probability"] > 0.00]
# Calculate the number of extant solution samples per
# category as a function of the reproduction probability
# and overall population size.
population_sizes = \
np.around(
(1.0 + survival_probability - natural_death_factor) * result_data["count"])
# Ensure that no population size drops below zero.
population_sizes = \
np.maximum(0.0, population_sizes)
# Augment the data frame so we can continue to keep
# track of the current population until we are done with
# this iteration for sure.
result_data = result_data.assign(
population_size=population_sizes)
# Convert the data types in the data frame to their most
# appropriate form.
result_data = result_data.convert_dtypes()
# Generate population for the next iteration.
#
# All I'm doing here in these two lines is just creating
# variables with shorter names. I'm not too familiar
# with pandas yet, so this conversion to a regular list
# could probably be done a lot less... badly.
solutions = result_data["solution"].array.to_numpy().tolist()
population_sizes = result_data["population_size"].array.to_numpy().tolist()
# Initialize the population variable.
population = [ ]
# Iterate over the solution and population size to be
# able to create the appropriate number of child
# elements for the next generation.
for solution, size in zip(solutions, population_sizes):
# Convert the population size of the given solution
# to an actual integer.
size = int(size)
# Generate the next generation's population.
#
# TODO: Figure out how to save population as a
# generator expression so we don't need to create
# the entire list in memory right away.
population += [
Solution(items, limit, solution) for _ in range(size)
]
# Create pairs of matched up population elements which
# will potentially undergo a crossover event.
pairings = combinations(population, 2)
# Initialize the pairings object from the generator.
pairings = list(pairings)
for A, B in pairings:
# Python strings are immutable, so we represent the
# temporary genome sequences as arrays, which later
# we'll "".join(a) and "".join(b) to come up with
# the actual strings.
a = [c for c in str(A)]
b = [d for d in str(B)]
# Simulate genetic mutations within population
# element A's genome.
for index, mutation in enumerate(rng.binomial(1, genetic_mutation_probability, len(items))):
if mutation:
a[index] = "0" if a[index] == "1" else "1"
# Simulate genetic mutations within population
# element B's genome.
for index, mutation in enumerate(rng.binomial(1, genetic_mutation_probability, len(items))):
if mutation:
b[index] = "0" if b[index] == "1" else "1"
# Make sure that a and b are not the same, since
# otherwise there is literally no point in doing a
# crossover event simulation.
if "".join(a) == "".join(b):
# Since both elements are made up of the same
# exact genetic components, we might as well
# skip this pair. After all, a crossover event
# between two identical elements would by
# definition be completely unnoticeable. Simply
# move on.
continue
# Now that we know that the two genetic elements are
# not identical, we can go ahead and check whether a
# cross-over even needs to be carried out. If it
# does, we go ahead and do that.
for c, d in zip(str(A), str(B)):
# We model the chance that a cross-over event
# will occur as a Bernoulli random variable with
# a probability of 0.57.
#
# Note, however, that this is the probability of
# a single cross-over event occuring. Multiple
# cross-over events could occur between two
# members of a population, albeit with
# significantly lower probability. In this model
# we consider cross-over events independent, and
# thus the probability of multiple cross-over
# events is equal to (0.57)^n, where n is equal
# to the number of cross-over events.
#
# Whether the assumption that cross-over events
# can actually be reasonably assumed to be
# independent is completely unbeknownst to me.
for index, crossover in enumerate(rng.binomial(1, crossover_probability, len(items))):
# Perform cross-over event on the current
# monomer, which we are metaphorically
# representing as a single character that
# could either be a zero or a one.
if crossover:
temp = a[index]
a[index] = b[index]
b[index] = temp
# Replace the solution values with the modified
# genome sequences.
A = Solution(items, limit, "".join(a))
B = Solution(items, limit, "".join(b))
# Add the current counts to the pandas series array.
for key, val in counter.items():
series[key] = series[key].append(Series(val))
# Initialize configuration variables.
maximum_x_value = None
minimum_x_value = 0
maximum_y_value = 1
minimum_y_value = 0
# Initialize the figure.
figure, axis = plt.subplots()
figure.set_dpi(300.0)
figure.set_size_inches(w=8.0, h=6.2)
# Configure the axis.
axis.axhline(y=0, linewidth=1.5, color="black")
axis.axvline(x=0, linewidth=1.5, color="black")
axis.grid(True, which="major")
# Plot all of the simulation series for which we have data.
for key in series.keys():
# Configure the plot data for this series.
data = series.get(key).to_list()[1:]
# The maximum x value should be the number of iterations
# performed.
if maximum_x_value is None or maximum_x_value < len(data):
maximum_x_value = len(data)
# The maximum y value should be set dynamically based
# on the values in the plot data.
if np.max(data) > maximum_y_value:
maximum_y_value = np.max(data) * 1.03
# Don't bother graphing uninteresting genomes. Only
# graph the ones that actually made it somewhere, in
# this case at least higher than eight.
if np.max(data) < 8:
# Move on to the next series.
continue
# Plot the graph.
axis.plot(data, label=f"{key}")
# Wrap up last-minute configuration of the graph.
axis.set_xlim(left=0, right=maximum_x_value)
axis.set_ylim(bottom=0, top=maximum_y_value + 1.02)
axis.set_title("Genetic Algorithm Simulation")
axis.set_xlabel("Iterations")
axis.set_ylabel("Population Count")
axis.xaxis.set_major_locator(MultipleLocator(1))
axis.xaxis.set_major_formatter("{x:.0f}")
axis.yaxis.set_major_locator(MultipleLocator(10))
axis.legend()
# Save the image just in case we want to keep it.
plt.savefig("simulation.png")
# Display the image using PyPlot.
plt.show()
