写在前面的话:
魔术方块就是中国所谓的九宫图(洛书):
他的特点就是横竖以及左右对角线的数字加起来都是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,以及方格数量。