python 遗传算法_遗传算法Python实战 008.魔术方块

写在前面的话：

魔术方块就是中国所谓的九宫图(洛书)：

他的特点就是横竖以及左右对角线的数字加起来都是15

一般来说，只要中间是5,其他的排列都可以通过头角互换(旋转和镜像)，比如：

6 7 2 
1 5 9
8 3 4

九宫格非常简单，就算是迭代，也能很快算出结果，如果要扩展成跟多的魔术方块，比如44或者干脆是1010呢？

？　？　？　？
？　？　？　？
？　？　？　？
？　？　？　？

我们通过算法来实现一下看看那

算法实现

首先导入需要的包，本次要用到一个新的方法：

import random,datetime
from bisect import bisect_left
from math import exp

bisect_left是Python内置的一个插入预测方法，他可以找到你需要插入数据所在的位置，功能特别好用，示例如下：

from bisect import bisect_left
a = list(range(10))
print(a)
# [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

#如果将3.6插入到列表a中，会出现在哪个位置。
idx = bisect_left(a,3.6)
print(idx)
# 4

然后初始化一堆参数：

• diagonalSize是你方块的数量，3就表示3*3的方块

• maxAge是本节算法需要着重说明的地方，后面解释

• nSquared 这个方块表面总共有多少个格子

• geneset是你的基因库，从1开始，到格子数量。

• expectedSum是对角线之和的值

• geneIndexes是每个基因的索引，与基因库里面的基因不同的是，索引从0开始

diagonalSize = 3
maxAge = 50
nSquared = diagonalSize * diagonalSize
geneset = [i for i in range(1, nSquared + 1)]
expectedSum = diagonalSize * (nSquared + 1) / 2
geneIndexes = [i for i in range(0, len(geneset))]

定义一个健壮性的类

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)

下面两个方法是用于健壮性检验的，其中：

• get_sums是内部方法，在get_fitness里面被调用，主要的作用是计算方块行列对角线的值分别是多少。

• get_fitness是健壮性建议函数，它把行列对角线所有数值都累加起来，然后计算与预期值的差距，如果等于0则表示完全符合预期(比如3*3的格子里面，行列对角线都是15)

def get_sums(genes):
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, southeastDiagonalSum

def get_fitness(genes):
rows, columns, northeastDiagonalSum, southeastDiagonalSum = \
        get_sums(genes)

    sumOfDifferences = sum(int(abs(s - expectedSum))
for s in rows + columns +
                           [southeastDiagonalSum, northeastDiagonalSum]
if s != expectedSum)

return Fitness(sumOfDifferences)

显示方法：解释略，除了显示当前方块的情况以外，还能够打印出行列对角线的和。

def display(candidate, startTime):
    timeDiff = datetime.datetime.now() - startTime

    rows, columns, northeastDiagonalSum, southeastDiagonalSum = \
        get_sums(candidate.Genes)

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)

进化方法：本节的进化方法又回到了最简单的状态，直接随机交换两个基因，完成进化

def mutate(chm):
    genes = chm.Genes[:]
    indexA, indexB = random.sample(geneIndexes,2)
    genes[indexA], genes[indexB] = genes[indexB], genes[indexA]
    return Chromosome(genes, get_fitness(genes))

种群类，多了一个成员变量Age——它的用途稍后解释

class Chromosome:
def __init__(self, genes, fitness):
self.Genes = genes
self.Fitness = fitness
self.Age = 0

本节算法的核心方法，进化过程这也是到现在为止，最复杂的进化方法，见代码中的注释：

#初始化初代种群的方法
def fnCustomCreate():
    gens = random.sample(geneset, len(geneset))
    ch1 = Chromosome(gens, get_fitness(gens))
return ch1

# 这是一个内部方法，用来不断生成进化种群
# 里面使用是Python特有的yield方法，用于返回生成的种群
# yield的用法，见Python语法原理。
def _get_improvement():
# 随机初始化初代种群，对于遗传算法来说，初始情况完全不重要
    parent = bestParent = fnCustomCreate()
yield bestParent
# 历史健壮性集合，这个集合可以记录历史上一些最佳的健壮性
    historicalFitnesses = [bestParent.Fitness]
while True:
        child = mutate(parent)
if parent.Fitness > child.Fitness:
# 本节算法用maxAge这个参数来控制单次迭代的数量
# 主要为了解决局部过大或者过小的问题，也就是说
if maxAge is None:
continue
            parent.Age += 1
# 初始化的时候设置maxAge为50，表示跟踪这个父本基因进化50次
# 生成出来的子代基因的健壮性还不如父本的话，这父本基因就直接被设置死亡
# 如果超过了代数，设置这个基因死亡，然后考虑换一个基因来进行替代。
# 划重点：这是模拟退火算法里面的核心思维
if maxAge > parent.Age:
continue
# 我们来计算，本次进化的健壮性，在整个历史进化过程中所有健壮性的位置。
            index = bisect_left(historicalFitnesses, child.Fitness, 0,
                                len(historicalFitnesses))
# 直接把这个位置转换成比例
            proportionSimilar = index / len(historicalFitnesses)
# 如果子代的健壮性很好，那么它的指数就会很高(如果排名历史第一(值为0)
# 则指数就会很高，百分之百会让这个子代变成下一次进化的父本。
# 反之，排名月底，成为父本的概率也就越低
# 就继续用这个父本进行进化
if random.random() < exp(-proportionSimilar):
                parent = child
continue
# 否则，用初始化的父本来替换成需要进化的父本
# 也就是丢弃本轮进化，从新回退到上一次最佳父本的状态
            bestParent.Age = 0
            parent = bestParent
continue
# 如果子体的健壮性要比父本好，则子体替换父本，成为新的进化的父本
# 替换了之后，子体的年龄要继承父本的年龄，并且继续增长。
if not child.Fitness > parent.Fitness:
            child.Age = parent.Age + 1
            parent = child
continue
        child.Age = 0
        parent = child
# 如果本次进化的子体比最佳父本效果还好，那么用这次进化的子体去替换最好的父本
# 然后把上一次最好的父本加入到历史父本中。
if child.Fitness > bestParent.Fitness:
            bestParent = child
yield bestParent
            historicalFitnesses.append(bestParent.Fitness)

def get_best():
    startTime = datetime.datetime.now()
    optimalFitness = Fitness(0)
# 迭代执行
# 直到到达进化目标
for improvement in _get_improvement():
        display(improvement,startTime)
if not optimalFitness > improvement.Fitness:
return improvement

模拟退火算法

退火是一种用于降低金属等材料内应力的方法。它的工作原理是通过加热，将金属升温到高温，然后让它慢慢冷却。我们都知道，热量会使金属膨胀。所以当金属膨胀的时候，金属材料之间的粘结会被松开，金属的分子结构就可以进行移动。这时，停止加热的时候，金属会缓慢冷却。当金属冷却后，金属的键再次收紧，杂质会被挤压出来，而且，直到他们相互之间会紧紧的粘合在一起，使之没有回旋的余地，从而减少了系统中的整体应力。

模拟退火是用遗传算法来突破局部最小或最大的过程。如果子代基因序列离目前的最佳解决方案很远，那么我们给它一个很高的概率继续进化下去。否则，那我们就做别的事情

执行结果如下：

# 最后一行是列的和，左边和右边是对角线的和
# 最后那个23，是健壮性
      [1, 4, 8] = 13
	  [3, 9, 6] = 18
	  [7, 2, 5] = 14
24 	  [11, 15, 19] 	 15
 - - - - - - - - - - - 23 0:00:00
 #……中间结果略
      [2, 9, 4] = 15
	  [7, 5, 3] = 15
	  [6, 1, 8] = 15
15 	 [15, 15, 15] 	 15
 - - - - - - - - - - - 0 0:00:05.995768
 # 0表示符合预期进化目标，完成进化

下面我们来修改一下我们的方块，比如改成4*4，修改基本参数如下：

diagonalSize = 4

执行结果如下:

      [5, 13, 3, 2] = 23
	  [7, 9, 14, 15] = 45
	  [6, 10, 4, 12] = 32
	  [16, 8, 11, 1] = 36
42 	 [34, 40, 32, 30] 	 19
 - - - - - - - - - - - 61 0:00:00
 # ……中间结果略

  [13, 4, 5, 12] = 34
  [16, 7, 10, 1] = 34
  [2, 9, 8, 15] = 34
  [3, 14, 11, 6] = 34
34 	 [34, 34, 34, 34] 	 34
 - - - - - - - - - - - 0 0:00:00.233368

最后来一个10*10的，修改基本参数如下：设置maxAge为5000

diagonalSize = 10
maxAge = 5000

这里运行效率的时候看人品，从几分钟到十几分钟不等(我这里这次用了2分多钟)，大家可以自行运行测试，运行结果如下：

      [4, 64, 10, 93, 19, 81, 60, 47, 59, 90] = 527
	  [16, 15, 24, 14, 21, 71, 48, 50, 40, 2] = 301
	  [70, 49, 34, 25, 23, 32, 57, 42, 37, 51] = 420
	  [63, 66, 27, 45, 13, 46, 54, 98, 82, 85] = 579
	  [7, 99, 6, 97, 87, 75, 33, 17, 20, 96] = 537
	  [86, 95, 61, 18, 77, 35, 56, 92, 94, 11] = 625
	  [26, 12, 79, 22, 1, 65, 36, 5, 29, 76] = 351
	  [67, 53, 100, 83, 39, 73, 62, 38, 89, 55] = 659
	  [41, 44, 91, 68, 84, 58, 3, 78, 72, 88] = 627
	  [52, 28, 9, 80, 30, 69, 31, 43, 74, 8] = 424
596 [432, 525, 441, 545, 394, 605, 440, 510, 596, 562] 	 374
 - - - - - - - - - - - 1896 0:00:00

 # 中间结果略……
 # 最后结果如下：行列对角线，全部等于505

      [96, 50, 25, 17, 28, 72, 61, 5, 75, 76] = 505
	  [66, 93, 18, 39, 83, 31, 78, 49, 8, 40] = 505
	  [34, 89, 52, 63, 62, 23, 64, 26, 24, 68] = 505
	  [16, 14, 27, 43, 60, 45, 32, 99, 74, 95] = 505
	  [47, 20, 56, 87, 88, 73, 19, 69, 42, 4] = 505
	  [86, 79, 12, 3, 77, 36, 54, 90, 58, 10] = 505
	  [22, 94, 97, 11, 1, 57, 46, 59, 48, 70] = 505
	  [55, 2, 98, 85, 44, 71, 33, 15, 67, 35] = 505
	  [30, 51, 82, 65, 21, 6, 37, 84, 29, 100] = 505
	  [53, 13, 38, 92, 41, 91, 81, 9, 80, 7] = 505
505  [505, 505, 505, 505, 505, 505, 505, 505, 505, 505] 	 505
 - - - - - - - - - - - 0 0:02:22.627786

有兴趣的同学，还可以尝试修改各种参数，比如maxAge，以及方格数量。