Magic squares are squares of numbers where each of the rows and columns and both major diagonals sum to the same value, and all of the numbers between 1 and n2 are used only once.
6 7 2
1 5 9
8 3 4We’ll start with the test stub.
import unittest
import datetime
import genetic
class MagicSquareTests(unittest.TestCase):
def test_size_3(self):
self.generate(3)In the test harness we’ll determine the range of numbers to use to fill the magic square from its diagonal size then calculate the expected sum. We then pass those to the fitness function.
def generate(self, diagonalSize):
nSquared = diagonalSize * diagonalSize
geneset = [i for i in range(1, nSquared + 1)]
expectedSum = diagonalSize * (nSquared + 1) / 2
def fnGetFitness(genes):
return get_fitness(genes, diagonalSize, expectedSum)
...For fitness we can count the number of rows, columns, and diagonals where the sum is equal to expectedSum.
def get_fitness(genes, diagonalSize, expectedSum):
rowSums = [0 for _ in range(diagonalSize)]
columnSums = [0 for _ in range(diagonalSize)]
southeastDiagonalSum = 0
northeastDiagonalSum = 0
for row in range(diagonalSize):
for column in range(diagonalSize):
value = genes[row * diagonalSize + column]
rowSums[row] += value
columnSums[column] += value
southeastDiagonalSum += genes[row * diagonalSize + row]
northeastDiagonalSum += genes[row * diagonalSize +
(diagonalSize - 1 - row)]
fitness = sum(1 for s in rowSums + columnSums +
[southeastDiagonalSum, northeastDiagonalSum]
if s == expectedSum)
return fitnessUnfortunately, that means if we want to display the sums we have to recalculate them. So, let’s extract a reusable function.
def get_sums(genes, diagonalSize):
rows = [0 for _ in range(diagonalSize)]
columns = [0 for _ in range(diagonalSize)]
southeastDiagonalSum = 0
northeastDiagonalSum = 0
for row in range(diagonalSize):
for column in range(diagonalSize):
value = genes[row * diagonalSize + column]
rows[row] += value
columns[column] += value
southeastDiagonalSum += genes[row * diagonalSize + row]
northeastDiagonalSum += genes[row * diagonalSize +
(diagonalSize - 1 - row)]
return rows, columns, northeastDiagonalSum, southeastDiagonalSumThen call it from the fitness function:
def get_fitness(genes, diagonalSize, expectedSum):
rows, columns, northeastDiagonalSum, southeastDiagonalSum = \
get_sums(genes, diagonalSize)
fitness = sum(1 for s in rows + columns +
[southeastDiagonalSum, northeastDiagonalSum]
if s == expectedSum)
return fitnessWe can call that function from display too.
def display(candidate, diagonalSize, startTime):
timeDiff = datetime.datetime.now() - startTime
rows, columns, northeastDiagonalSum, southeastDiagonalSum = \
get_sums(candidate.Genes, diagonalSize)
for rowNumber in range(diagonalSize):
row = candidate.Genes[
rowNumber * diagonalSize:(rowNumber + 1) * diagonalSize]
print("\t ", row, "=", rows[rowNumber])
print(northeastDiagonalSum, "\t", columns, "\t", southeastDiagonalSum)
print(" - - - - - - - - - - -", candidate.Fitness, timeDiff)This will produce output like the following:
[5, 1, 8] = 14
[6, 4, 9] = 19
[3, 7, 2] = 12
15 [14, 12, 19] 11
- - - - - - - - - - - 1 0:00:00.001000Now, since one of the constraints is that we have to use all the numbers between 1 and n2, we’ll pass a custom_create function that produces a random permutation of all values in that range
import random
...
def generate(self, diagonalSize):
...
def fnCustomCreate():
return random.sample(geneset, len(geneset))
optimalValue = 2 + 2 * diagonalSize
startTime = datetime.datetime.now()
best = genetic.get_best(fnGetFitness, nSquared, optimalValue,
geneset, fnDisplay, fnMutate,
fnCustomCreate)We’ll also pass a companion function for custom_mutate so we don’t have to worry about duplication being introduced by the engine’s built-in mutate. We can simply swap two genes but for performance we’re going to pass in the index options. This can often provide a good performance boost versus generating a list of indexes each time we enter the function.
def mutate(genes, indexes):
indexA, indexB = random.sample(indexes, 2)
genes[indexA], genes[indexB] = genes[indexB], genes[indexA]Here’s the full test harness.
def generate(self, diagonalSize):
nSquared = diagonalSize * diagonalSize
geneset = [i for i in range(1, nSquared + 1)]
expectedSum = diagonalSize * (nSquared + 1) / 2
def fnGetFitness(genes):
return get_fitness(genes, diagonalSize, expectedSum)
def fnDisplay(candidate):
display(candidate, diagonalSize, startTime) geneIndexes = [i for i in range(0, len(geneset))]
def fnMutate(genes):
mutate(genes, geneIndexes)
def fnCustomCreate():
return random.sample(geneset, len(geneset))
optimalValue = 2 + 2 * diagonalSize
startTime = datetime.datetime.now()
best = genetic.get_best(fnGetFitness, nSquared, optimalValue,
geneset, fnDisplay, fnMutate,
fnCustomCreate)
self.assertTrue(not optimalValue > best.Fitness)And, when we run it we get a result like the following:
...
[8, 3, 4] = 15
[1, 9, 5] = 15
[2, 7, 6] = 15
15 [11, 19, 15] 23
- - - - - - - - - - - 5 0:00:00.001001We can run the test many times but it will rarely find a valid solution.
The problem is the fitness function is written exclusively rather than inclusively. It only gives credit when the sum exactly matches the expected sum, so there’s no partial credit for the genetic engine to take advantage of. This is a problem we’ve encountered before. The solution is to take the sum of the differences between the actual sum and the expected sum for each row, column, and diagonal. That makes zero optimal, so we’ll need a Fitness object to help reverse the greater-than logic used to compare fitnesses in the engine.
class Fitness:
def __init__(self, sumOfDifferences):
self.SumOfDifferences = sumOfDifferences
def __gt__(self, other):
return self.SumOfDifferences < other.SumOfDifferences
def __str__(self):
return "{}".format(self.SumOfDifferences)Next we update get_fitness:
def get_fitness(genes, diagonalSize, expectedSum):
rows, columns, northeastDiagonalSum, southeastDiagonalSum = \
get_sums(genes, diagonalSize)
sumOfDifferences = sum(int(abs(s - expectedSum))
for s in rows + columns +
[southeastDiagonalSum, northeastDiagonalSum]
if s != expectedSum)
return Fitness(sumOfDifferences)And the optimal fitness value in our test.
def generate(self, diagonalSize):
...
optimalValue = Fitness(0)
startTime = datetime.datetime.now()Now when we run the test we get a valid solution about 60% of the time but only because it just happens to find a sequence of swaps that drive improvement. When it doesn’t find a sequence
[4, 8, 3] = 15
[2, 6, 7] = 15
[9, 1, 5] = 15
18 [15, 15, 15] 15
- - - - - - - - - - - 3 0:00:00.002005it gets stuck at a position where it requires at least 2 swaps to make progress. We’re hitting a local minimum.
To fix the local minimum/local maximum issue we’re going to allow the current genetic line to die out. The first step in that is tracking how many generations have passed since the last improvement. We’ll call this its Age.
class Chromosome:
def __init__(self, genes, fitness):
self.Genes = genes
self.Fitness = fitness
self.Age = 0Next, we’ll add an optional parameter that allows us to set an upper limit on the age of a genetic line. That parameter will be passed through to _get_improvement.
def get_best(get_fitness, targetLen, optimalFitness, geneSet, display,
custom_mutate=None, custom_create=None, maxAge=None):
...
for improvement in _get_improvement(fnMutate, fnGenerateParent, maxAge):For the next part we need to import a couple of module functions.
from bisect import bisect_left
from math import expThen we need to separate the best parent from the current parent so that we still have something to compare to when a genetic line dies out. We’re also going to keep a list of the fitnesses of the historical best parents.
def _get_improvement(new_child, generate_parent, maxAge):
parent = bestParent = generate_parent()
yield bestParent
historicalFitnesses = [bestParent.Fitness]
while True:Next, we want to make sure we retain the current functionality if maxAge is not provided.
while True:
child = new_child(parent)
if parent.Fitness > child.Fitness:
if maxAge is None:
continueHowever, when the child’s fitness is worse than that of its parent, the most traveled path through the code, and maxAge is provided, then we need to check whether or not the genetic line’s age has reached the maximum.
continue
parent.Age += 1
if maxAge > parent.Age:
continueIf so, we may allow the genetic line to die and replace it with something else. We’re going to do that with simulated annealing.
We need to determine how far away the child’s fitness is from the best fitness. If we had a numeric fitness it would be a simple calculation. But our fitnesses aren’t always numeric. No problem, we’ll figure out where it would be in the realm of historical fitnesses.
index = bisect_left(historicalFitnesses, child.Fitness, 0, len(historicalFitnesses))Then get its proximity to the best fitness.
difference = len(historicalFitnesses) - indexWhether calculating the difference directly from the fitness, or by using its index into the historical fitnesses as we do here, we then convert it to a proportion by dividing it by the best value, in this case the highest index.
proportionSimilar = difference / len(historicalFitnesses)Then we raise Euler’s number, roughly 2.718, to the power of the proportion negated. The result is a floating point number that approaches 1 if the child’s fitness is far away from the best fitness, but approaches ~0.36 when its fitness is very close to the best fitness. In the following sample chart assume the maximum current index is 50:
index | difference | proportion similar | e^-proportion ------------------------------------------------------- 0 | 50 | 0.0 | 1.00 5 | 45 | 0.1 | 0.90 10 | 40 | 0.2 | 0.82 40 | 10 | 0.8 | 0.45 45 | 5 | 0.9 | 0.41 50 | 0 | 1.0 | 0.37
A child whose fitness is close to the current best will have a high index (because the higher/better fitnesses are at the end of the array) and a low difference from the best fitness. As a result, it will have a lower chance of becoming the new parent. A child that has a fitness far from the current best will have a low index and high difference, and thus a high chance of becoming the new parent.
Next we pick a random number and if that random number is smaller than e-proportion then the child becomes the new parent.
if random.random() < exp(-proportionSimilar):
parent = child
continueOtherwise we replace the parent with the best parent and reset its age to zero so it has time to anneal.
continue
parent = bestParent
parent.Age = 0
continueNext we consider what to do if the child’s fitness is not lower than that of its parent.
continue
if not child.Fitness > parent.Fitness:
# same fitness
child.Age = parent.Age + 1
parent = child
continueWhen the child has a better fitness than its parent, we reset its age to zero and make it the new parent.
continue
parent = child
parent.Age = 0Finally, when we find a child whose fitness is better than that of the best parent, we replace the best parent and append its fitness to the list of historical fitnesses.
parent.Age = 0
if child.Fitness > bestParent.Fitness:
yield child
bestParent = child
historicalFitnesses.append(child.Fitness)Now, to use it in our test harness we just need to set maxAge. But what value should we use? How do we pick a good maximum age? Let’s sample the average run times for the projects where we’ve encountered local minimums or maximums.
| 50 | 500 | 5000 | no max ----------------------------------------------------------------------------- Knights | 0.63 +/- 0.47 | 0.68 +/- 0.52 | 0.66 +/- 0.53 | 0.61 +/- 0.46 Magic Square | 0.01 +/- 0.01 | 0.04 +/- 0.06 | 0.39 +/- 0.47 | *
That’s interesting. The Knights project doesn’t appear to benefit from allowing genetic lines to age out. Magic squares, however, clearly benefits not only from having a maximum age but also from having a low one.
def test_size_3(self):
self.generate(3, 50) def generate(self, diagonalSize, maxAge):
...
best = genetic.get_best(fnGetFitness, nSquared, optimalValue,
geneset, fnDisplay, fnMutate,
fnCustomCreate, maxAge)Now the test can quickly find a solution for a size 3 magic square every time.
[6, 7, 2] = 15
[1, 5, 9] = 15
[8, 3, 4] = 15
15 [15, 15, 15] 15
- - - - - - - - - - - 0 0:00:00.008000Size-5 magic squares need a higher maximum age.
def test_size_5(self):
self.generate(5, 500)It can find a magic square every time but we’re really giving the simulated annealing code a workout to get it.
[25, 3, 10, 4, 23] = 65
[9, 21, 8, 14, 13] = 65
[22, 6, 5, 15, 17] = 65
[2, 16, 24, 12, 11] = 65
[7, 19, 18, 20, 1] = 65
65 [65, 65, 65, 65, 65] 64
- - - - - - - - - - - 1 0:00:00.285760
[13, 16, 7, 4, 25] = 65
[22, 19, 5, 9, 10] = 65
[20, 3, 17, 14, 11] = 65
[2, 6, 24, 15, 18] = 65
[8, 21, 12, 23, 1] = 65
65 [65, 65, 65, 65, 65] 65
- - - - - - - - - - - 0 0:00:00.921482 def test_size_10(self):
self.generate(10, 5000)Given time, it can build larger squares too:
[31, 9, 50, 21, 81, 92, 27, 34, 97, 63] = 505
[ 8, 80, 86, 29, 40, 24, 100, 55, 72, 11] = 505
[87, 53, 69, 95, 15, 67, 94, 6, 17, 2] = 505
[79, 30, 37, 28, 38, 71, 42, 96, 62, 22] = 505
[54, 46, 56, 93, 3, 74, 59, 52, 23, 45] = 505
[35, 89, 98, 83, 32, 4, 10, 5, 58, 91] = 505
[19, 75, 1, 41, 61, 70, 33, 76, 47, 82] = 505
[20, 44, 78, 12, 85, 36, 65, 99, 18, 48] = 505
[88, 13, 14, 64, 73, 60, 49, 25, 68, 51] = 505
[84, 66, 16, 39, 77, 7, 26, 57, 43, 90] = 505
505 [505, 505, 505, 505, 505, 505, 505, 505, 505, 505] 505
- - - - - - - - - - - 0 0:04:24.477191We saw that using various maximum ages on the Knight’s project had no particular effect. What do you think will happen if you use a maximum age for the other projects? Try it.
We’ll benchmark this project with size-4 magic squares because size-5 runs a bit slower than I like in a benchmark.
def test_size_4(self):
self.generate(4, 50)
def test_benchmark(self):
genetic.Benchmark.run(self.test_size_4)We rewrote a core function of the genetic module so we’ll update the benchmarks to make sure there was no unintended slowdown in the ability to solve the previous projects.
Updated Benchmarks
project | average (seconds) | standard deviation
-----------------------------------------------------------
Guess Password | 1.31 | 0.41
One Max | 1.22 | 0.17
Sorted Numbers | 1.23 | 0.83
Queens | 1.39 | 1.03
Graph Coloring | 0.78 | 0.34
Cards | 0.01 | 0.01
Knights | 0.61 | 0.46
Magic Square | 0.33 | 0.56
In this chapter we encountered the first project where our forward progress was significantly impacted by a local minimum. We resolved this by adding an option to the engine that allows it to escape local minimums/maximums through the use of simulated annealing. This is also the first time we passed the gene indexes to the mutate function. We’ll be using this again.
The final code for this chapter is available from:
https://github.com/handcraftsman/GeneticAlgorithmsWithPython/tree/master/ch08