遺傳算法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)
下面兩個方法是用於健壯性檢驗的，其中：
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，以及方格數量。