1. 遺傳程式設計簡介

0x1：什麼是遺傳程式設計演算法，和傳統機器學習演算法有什麼區別

傳統上，我們接觸的機器學習演算法，都是被設計為解決某一個某一類問題的確定性演算法。對於這些機器學習演算法來說，唯一的靈活性體現在引數搜尋空間上，向演算法輸入樣本，演算法藉助不同的優化手段，對引數進行調整，以此來得到一個對訓練樣本和測試樣本的最佳適配引數組。

遺傳程式設計演算法完全走了另一外一條路，遺傳程式設計演算法的目標是編寫一個程度，這個程式會嘗試自動構造出解決某一問題的最佳程度。從本質上看，遺傳程式設計演算法構造的是一個能夠構造演算法的演算法。

另一方面，我們曾經討論過遺傳演算法，遺傳演算法是一種優化技術，就優化技術而言，無論是何種形式的優化，演算法或度量都是預先設定好的，而優化演算法所做的工作就是嘗試為其找到最佳引數。和優化演算法一樣，遺傳程式設計也需要一種方法來度量題解的優劣程度。但與優化演算法不同的是，遺傳程式設計中的題解並不僅僅是一組用於給定演算法的引數，相反，在遺傳程式設計中，連同演算法本身及其所有引數，都是需要搜尋確定的。

從某種程度上來說，遺傳程式設計和遺傳演算法的區別在於，進化的基本單位不同，

• 遺傳優化：進化的基本單位是模型可變引數
• 遺傳程式設計：進化的基本單位是新演算法以及新演算法的引數

0x2：遺傳程式設計和進化論的關係

遺傳演算法是受達爾文的進化論的啟發，借鑑生物進化過程而提出的一種啟發式搜尋演算法，因此遺傳演算法 ( GA , Genetic Algorithm ) 也稱進化演算法 。 因此，在討論遺傳程式設計的時候，會大量借用進化論中的術語和概念，為了更好地討論遺傳演算法，我們先介紹一些基本生物進化概念，

• 基因 ( Gene )：一個遺傳因子，種群中的最基本單元。
• 染色體 ( Chromosome )：一組的基因。
• 個體 ( individual )：單個生物。在遺傳演算法中，個體一般只包含一條染色體。
• 種群 ( Population )：由個體組成的群體。生物的進化以種群的形式進化。
• 適者生存 ( The survival of the fittest )：對環境適應度高的個體參與繁殖的機會比較多，後代就會越來越多。適應度低的個體參與繁殖的機會比較少，後代就會越來越少。

生物所處的環境起到一個提供生存壓的作用（反饋），雖然縱觀整個地球歷史，環境的因素是在不斷變化的（有時甚至變化的還很快），但是在某個時間段內（例如5000年內）是基本保持不變的，而物種進化的目的就是通過一代代的繁衍，逐漸適應（擬合）當前的環境，並和其他物種達到最優平衡（納什均衡）。

遺傳程式設計演算法就是模擬了生物進化的過程，簡單說來說，

• 生物進化的環境由一個使用者定義的任務（user-defined task）所決定，演算法由一組初始的題解（程式）開始展開競爭。這裡所謂的任務可以是多種形式，
  • 一種競賽（game）：各個題解（程式）在競賽中直接展開競爭
  • 個體測試：測出哪個題解（程式）的執行效果更好
• 遺傳演算法將基因抽象為題解中最小的隨機變數因子（例如模型中的可變引數）
• 一個問題的解由很多這樣的隨機變化因子組成，演算法將問題的解編碼成個體的染色體（染色體是基因的集合）
• 單個個體包含若干個染色體，個體包含的染色體（題解）越多和越好，則個體的適應度就越好。在實際工程中，為了簡化演算法，常常假設一個個體只有一條染色體
• 多個個體組成種群，種群中適應度（Fitness）高的個體獲得較高概率的繁殖機會，從而導致適應度高的個體會越來越多，經過N代的自然選擇後，儲存下來的個體都是適應度很高的
• 繁殖過程中，演算法會評估並挑選出本輪表現最好的一部分題解題解（程式），並對程式的某些部分以隨機（一定概率）的方式進行修改，包括：　　　　
  • 基因交叉（Acrossover）：在最優題解之間，挑選部分隨機變數因子進行彼此互換。遺傳演算法交叉比人體內染色體交叉要簡單許多。遺傳演算法的染色體是單倍體，而人體內的真正的染色體是雙倍體。下圖是遺傳演算法中兩條染色體在中間進行交叉的示意圖，
  • 基因突變（Mutation）：在最優題解上，直接對某些隨機變數因子（基因位）進行隨機修改。下圖是遺傳演算法中一條染色體在第二位發生基因變異的示意圖，
• 經過繁殖過程，新的種群（即新的一組解）產生，稱為“下一代”，理論上，這些新的題解基於原來的最優程式，但又不同於它們。這些新產生的題解和舊的最優題解會一起進入下一輪自然選擇階段
• 上述繁殖過程重複多次，直到達到收斂條件，包括，
  • 找到了全域性最優解
  • 找到了表現足夠好的解
  • 題解在歷經數代之後都沒有得到任何改善
  • 繁衍的代數達到了規定的限制
• 最終，歷史上適應度最高個體所包含的解，作為遺傳演算法的輸出

下圖是遺傳演算法的流程圖，

0x3：遺傳程式設計的不同型別

從大的方面看，遺傳程式設計的兩個重要概念是基因型和表現型，

• 基因型就是種群個體的編碼；
• 表現型是種群個體所表示的程式片段；

其實遺傳演算法領域的研究中，這兩個方面的研究都有，但是，因為遺傳程式設計很難直接處理程式片段（表現型）（例如一段C++可執行程式碼、或者是一段python程式碼），因為基於隨機變異得到的新程式碼很可能無法通過編譯器語法檢查。

但是相比之下，遺傳演算法反而容易處理程式片段的內在結構（基因型）（例如C++程式碼的AST抽象語法樹）。

所以，筆者認為基因型的遺傳演算法研究才是更有研究價值的一個方向，本文的討論也會圍繞基因型的遺傳演算法展開。

根據基因型形態的不同，遺傳程式設計方法可以分為三種：

• 線性遺傳程式設計
• 基於樹的遺傳程式設計
• 基於圖的遺傳程式設計

1. 線性遺傳程式設計

線性遺傳程式設計有廣義和狹義之分，

• 廣義線性遺傳程式設計將候選程式編碼進定長或者變長的字串，即基因型是線性字串，包括
  • Multi-Expression Programming (MEP)
  • Grammatical Evolution (GE)
  • Gene Expression Programming (GEP)
  • Cartesian Genetic Programming (CGP)：該演算法是一種很適合電路設計的遺傳程式設計演算法，比如我們要用兩個加操作兩個減操作和兩個乘操作得到如下運算，
    • 笛卡爾遺傳程式設計將下面的一個候選程式編寫進字串"001 100 131 201 044 254 2573"。字串中的三位數字“xyz"表示x操作的輸入是y和z兩個連線，字串中最後的四位數字"opqr"表示輸出opqr四個連線。笛卡爾遺傳程式設計只用變異操作，而不用交叉操作。
  • Genetic Algorithm for Deriving Software (GADS)
• 狹義線性遺傳程式設計中的候選程式是組合語言或者高階程式語言程式（例如C程式）。一個狹義線性遺傳程式設計的個體可以是一段簡單 C 語言指令，這些指令作用在一定數量預先定義的變數或者常量上（變數數量一般為指令個數的4倍）。下圖是一個狹義線性遺傳程式設計候選程式的示例，

• ，可以看到，變數數量和指令數量都是固定的，通過不同的排列組合方式得到不同的程式碼表現形式

http://www.doc88.com/p-630428999834.html
https://pdfs.semanticscholar.org/958b/f0936eda72c3fc03a09a0e6af16c072449a1.pdf

2. 基於樹的遺傳程式設計 

基於樹的遺傳程式設計的基因型是樹結構。基於樹的遺傳程式設計是遺傳程式設計最早的形態，也是遺傳程式設計的主流方法。

大多數程式語言，在編譯或解釋時，首先會被轉換成一棵解析樹（Lisp程式語言及其變體，本質上就是一種直接訪問解析樹的方法），例如下圖，

樹上的節點有可能是枝節點也可能是葉節點，

• 枝節點代表了應用於其子節點之上的某一種操作
• 葉節點代表了某個引數或常量值

例如上圖中，圓形節點代表了應用於兩個分支（Y變數和常量值5）之上的求和操作。一旦我們求出了此處的計算值，就會將計算結果賦予上方的節點處。相應的，這一計算過程會一直向下傳播，直到遍歷所有的葉子節點（深度優先遞迴遍歷）。 

如果對整棵樹進行遍歷，我們會發現它相當於下面這個python函式：

在遺傳變異方面，基於樹的遺傳程式設計的演化操作有兩種，變異和交叉，

• 變異：基於樹的遺傳程式設計的變異操作有兩種（區別在於變異的範圍不同），
  • 一種是隨機變換樹中的符號或者操作符
  • 另一種是隨機變換子樹
  • ，該圖左下角是變換符號或者操作符的結果，右下角是變換子樹的結果。 
• 交叉：兩個顆樹之間隨機交換子樹
  • ，兩棵樹之間的部分節點發生了隨機互換

3. 基於圖的遺傳程式設計

樹是一種特殊的圖，因此人們很自然地想到將基於樹的遺傳程式設計擴充套件到基於圖的遺傳程式設計。下圖就是基於圖的遺傳程式設計的基因型的一個示例。 

Relevant Link: 

《Adaptation in Natural and Artificial Systems》 John Henry Holland 1992
http://www.algorithmdog.com/%e9%81%97%e4%bc%a0%e7%ae%97%e6%b3%95%e7%b3%bb%e5%88%97%e4%b9%8b%e4%b8%80%e9%81%97%e4%bc%a0%e7%ae%97%e6%b3%95%e7%ae%80%e4%bb%8b
《Evolving Evolutionary Algorithms using Linear Genetic Programming (2005)》
《A comparison of several linear genetic programming techniques》Oltean, Mihai, and Crina Grosan.  Complex Systems 14.4 (2003): 285-314.
https://www.jianshu.com/p/a953066cb2eb 

2. 遺傳程式設計的數學基礎

這個章節，我們用數學的形式化視角，來重新審視一下遺傳演算法。

0x1：基本數學符號定義 

因為遺傳演算法中有各種各樣的編碼方式、變異操作、交叉操作和選擇操作，遺傳演算法的形態呈現多樣性。

為了簡化分析，我們這裡假設一個典型遺傳演算法，即，

• 編碼方式是二進位制編碼：基因的取值只能是0或者1
• 變異操作將所有染色體所有基因位以恆定 pm 的概率翻轉
• 交叉操作選擇選擇相鄰的個體，以 pc 的概率決定是否需要交叉。如果需要交叉，隨機選擇一個基因位，並交換這個基因位以及之後的所有基因
• 每一代的新種群選擇操作採用輪盤賭演算法（依據概率大小）：有放回地取樣出原種群大小的新一代種群，個體 I_i 的取樣概率如下所示，
  •  

0x2：模式定理 - 概率視角看基因模式的遺傳

模式定理是遺傳演算法創始人 J.Holland 在其突破性著作《Adaptation in Natural and Artificial Systems》引入的，用於分析遺傳演算法的工作原理。

模式是指基因編碼空間中，由一類相似的基因組抽象得到的pattern，比如 [0，*,*,1] 就是一個模式。染色體[0,1,0,1]和[0,0,0,1]都包含該模式。

在具體討論模式定理之前，我們先介紹一些符號，

【模式定理】

在本章定義的典型遺傳演算法中，下面公式成立：

這個公式看起來有點複雜，其實非常簡單，我們逐個部分來分解，

• 選擇操作對模式H在下一代出現的影響是固定的，即：

• 某個模式在繁衍中，既有可能發生變異，也有可能發生交叉，所以有後面兩個括號相乘
• 某個模式在變異中，變異操作將所有基因位以 pm 的概率翻轉，因此模式H不被破壞的概率為(1−pm)^O(H)。當0<=x<=1和n=1,...時，不等式(1−pm)^O(H) >= 1−O(H)∗pm成立，從而經過變異操作，模式H的出現概率為，
• 某個模式在交叉中，交叉操作選擇選擇相鄰的個體，以 pc 的概率決定是否需要交叉。如果需要交叉，隨機選擇一個基因位，並交換這個基因位以及之後的所有基因。因此模式H不被破壞的概率為(1−pc)(1−L(H)/L−1) >= 1 − pc∗L(H)/L−1。經過交叉操作，模式H的出現概率為，

總體來說，遺傳演算法需要在，選擇操作引起的收斂性和變異交叉操作引起的多樣性之間取得平衡。

模式定理的通俗說法是這樣的，低階、短序以及平均適應度高於種群平均適應度的模式在子代中呈指數增長。

低階、短長以及平均適應度高於種群平均適應度的模式H,

此時，

即模式H呈現指數增長。

0x3：馬爾柯夫鏈分析 - 遺傳程式設計收斂性分析

這個小節我們來討論一個有些理論化的問題，即：遺傳程式設計演算法經過一定次數的迭代後，是否會收斂到某個穩態？如果會達到穩態，遺傳程式設計的收斂速度是多少？

要解決這個問題，我們需要引入數學工具，馬爾柯夫鏈，有如下定義。

• 用 p^t 表示第 t 時刻的不同狀態的概率
• P 表示轉移概率矩陣，其中 P_i,j表示從第 i 個狀態轉移到第 j 個狀態的概率
• 齊次馬爾科夫鏈的第 t+1 時刻的狀態只和第 t 時刻有關，可以用公式 p^t+1=p^tP 表示
• 若存在一個自然數 k，使得 P^k 中的所有元素大於0，則稱 P 為素矩陣。隨著 k 趨近於無窮，P^k 收斂於 P^∞=1^Tp^∞, 其中p^∞=p⁰ lim_k→∞P^k=p⁰ 是和初始狀態無關的唯一值，並且所有元素大於0。這其實是由馬爾柯夫鏈穩態定理決定的。 

我們把整個種群的狀態看成馬爾科夫鏈的一個狀態 s，交叉、變異和選擇操作則構建了一個概率轉移矩陣。一般情況下，0<pm<1，0<=pc<=1，即物種變異一定會發生，但不是必然100%發生。我們來分析一下這時的概率轉移矩陣的性質。

• 讓 C，M，S 分別表示交叉、變異和選擇操作帶來的概率轉移，整體概率轉移矩陣 P=CMS
  • 經過變異操作，種群狀態 s_i 轉化成種群狀態 s_j 的概率 M_i,j=(pm)^h(1−pm)^n∗l-h>0，其中h是兩個種群之間不同值的基因位數量。也就是說，M 是素矩陣
  • 經過選擇操作，種群狀態 s_i 保持不變的概率，也就是說, S 的所有列必定有一元素大於0。我們也可以知道概率轉移矩陣 P 是素矩陣

標準的優化演算法分析第一個要關心的問題是，優化演算法能不能收斂到全域性最優點。假設全域性最優點的適應度值為maxf，收斂到全域性最優點的定義如下，

一言以蔽之，典型遺傳演算法並不收斂。

根據概率轉移矩陣收斂定理，我們可以知道典型遺傳演算法會收斂到一個所有種群狀態概率都大於0的概率分佈上（穩態）。因此之後，不包含全域性最優解的種群一定會不停出現，從而導致上面的公式不成立。

但是筆者這裡要強調的是，這章討論的典型遺傳演算法在實際工程中是幾乎不存在的，實際上，幾乎所有遺傳演算法程式碼都會將保持已發現最優解。加了這個變化之後的遺傳演算法是收斂的。

還是根據上述概率轉移矩陣收斂定理，我們可以知道遺傳演算法會收斂到一個所有種群狀態概率都大於0的概率分佈上，那麼包含全域性最優解的種群一定會不停出現，保持已發現最優解的做法會使得上面的公式成立。

Relevant Link: 

Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence
Rudolph, Günter. “Convergence analysis of canonical genetic algorithms.” Neural Networks, IEEE Transactions on 5.1 (1994): 96-101.
http://www.algorithmdog.com/%e9%81%97%e4%bc%a0%e7%ae%97%e6%b3%95%e7%b3%bb%e5%88%97%e4%b9%8b%e4%b8%89%e6%95%b0%e5%ad%a6%e6%91%86%e6%91%86%e6%89%8b%ef%bc%8c%e5%be%88%e6%83%ad%e6%84%a7%ef%bc%8c%e5%8f%aa%e5%81%9a%e4%ba%86

3. 典型遺傳演算法的一些變種衍生演算法

自 John Henry Holland 在1992年提出《Adaptation in Natural and Artificial Systems》論文後，遺傳程式設計又得到了大量研究者的關注和發展，提出了很多改進型的衍生演算法。雖然這些演算法在工業場景中不一定都適用，但是筆者覺得我們有必要了解和學習一下它們各自的演算法思想，有利於我們在遇到各自的實際問題的時候，舉一反三。

0x1：交叉變種

典型交叉變異隨機選擇兩條染色體，按照pc的概率決定是否交叉，如果選擇交叉則隨機選擇一點並將這點之後的基因交換。這種交叉方式被稱為單點雜交。

1. 多點雜交

多點雜交指定了多個交換點用於父本的基因交換重組，具體的執行過程如下圖所示，

多點雜交改進的是突變率。 

2. 均勻雜交

單點和多點雜交演算法存在一個問題，雜交的染色體中某些部分的基因會被過早地捨棄，這是由於在交換前它們必須確定交換父本染色體交換位前面還是後面的基因，從而對於那些無關的基因段在交換前就已經收斂了。

均勻雜交演算法（Uniform Crossover）就可以解決上述演算法的這種侷限性，該演算法的主要過程如下：

• 首先隨機選擇染色體上的交換位
• 然後隨機確定交換的基因是父本染色體上交換位的前部分基因，還是後部分基因（隨機過程）
• 最後對父本染色體的基因進行重組從而產生新的下一代個體

3. 洗牌雜交

洗牌雜交的最大特點是通常將染色體的中點作為基因的交換點，即從每個父本中取它們一半的基因重組成新的個體。

另外針對於實值編碼方式，還有離散雜交、中間雜交、線性雜交和擴充套件線性雜交等演算法。

0x2：選擇策略變種

精英保留策略是一種典型的選擇策略。精英保留策略是指每次迭代都保留已發現的最優解。這個策略是顯而易見的，我們不可能捨棄已發現的最優解，而只使用最後一代種群的最優解。同時，採用精英保留策略的典型遺傳演算法是保證收斂到全域性最優解的。

1. 輪盤賭選擇策略

輪盤賭選擇策略是基於概率進行選擇策略。輪盤賭演算法有放回地取樣出原種群大小的新一代種群，個體 I_i 的取樣概率如下所示，

從概率上看，在某一輪中，即使是適應度最差的個體，也存在一定的機率能進入到下一輪，這種策略提高了多樣性，但減緩了收斂性。

2. 錦標賽選擇策略

錦標賽法從大小為 n 的種群隨機選擇 k(k小於n) 個個體，然後在 k 個個體中選擇適應度最大的個體作為下一代種群的一個個體。反覆多次，直到下一代種群有 n 個個體。 

0x3：種群繁衍策略變種 - 多種群並行

在大自然，物種的進化是以多種群的形式併發進行的。一般來說，一個物種只有一個種群了，意味著這個物種有滅亡的危險（例如恐龍）。

受此啟發，人們提出了多種群遺傳演算法。多種群遺傳演算法保持多個種群同時進化，具體流程如下圖所示，

多種群遺傳演算法和遺傳演算法執行多次的區別在於移民，種群之間會通過移民的方式交換基因。這種移民操作會帶來更多的多樣性。

0x4：自適應遺傳演算法

遺傳演算法中，決定個體變異長度的主要因素有兩個：交叉概率pc，和變異概率pm。

在實際工程問題中，需要針對不同的優化問題和目標，反覆實驗來確定pc和pm，調參成本很高。

而且在遺傳演算法訓練的不同階段，我們需要不同的pc和pm，

• 當種群中各個個體適應度趨於一致或者趨於區域性最優時，使pc和pm增加，增加擾動。使得種群具有更大的多樣性，跳出潛在的區域性最優陷阱
• 當群體適應度比較分散時，使pc和pm減少。使得適應度高的個體和適應度低的個體保持分開，加快種群收斂速度
• 不同個體也應該有不同的pc和pm：
  • 對於適應度高的個體，我們應該減少pc和pm以保護他進入下一代
  • 反之對適應度低的個體，我們應該增加pc和pm以增加擾動，提高個體多樣性

Srinivas.M and Patnaik.L.M (1994) 為了讓遺傳演算法具備更好的自適應性，提出來自適應遺傳演算法。在論文中，pc和pm的計算公式如下：

0x5：混合遺傳演算法 

遺傳演算法的全域性搜尋能力強，但區域性搜尋能力較弱。這句話怎麼理解呢？

比如對於一條染色體，遺傳演算法並不會去看看這條染色體周圍區域性的染色體適應度怎麼樣，是否比這條染色體好。遺傳演算法會通過變異和交叉產生新的染色體，但新產生的染色體可能和舊染色差的很遠。因此遺傳演算法的區域性搜尋能力差。

相對的，梯度法、爬山法和貪心法等演算法的區域性搜尋能力強，運算效率也高。

受此啟發，人們提出了混合遺傳演算法，將遺傳演算法和這些演算法結合起來。混合遺傳演算法的框架是遺傳演算法的，只是生成新一代種群之後，對每個個體使用區域性搜尋演算法尋找個體周圍的區域性最優點。

總體來說，遺傳演算法和梯度法分別代表了隨機多樣性優化和漸進定向收斂性優化的兩種思潮，取各自的優點是一種非常好的思路。

Relevant Link: 

Srinivas M, Patnaik L M. Adaptive probabilities of crossover and mutation in genetic algorithms[J]. Systems, Man and Cybernetics, IEEE Transactions on, 1994, 24(4): 656-667. 
http://www.algorithmdog.com/%e9%81%97%e4%bc%a0%e7%ae%97%e6%b3%95%e7%b3%bb%e5%88%97%e4%b9%8b%e5%9b%9b%e9%81%97%e4%bc%a0%e7%ae%97%e6%b3%95%e7%9a%84%e5%8f%98%e7%a7%8d 

4. 用遺傳程式設計自動生成一個能夠擬合特定資料集的函式

0x1：用多項式迴歸擬合一個數據集

這個章節，我們來完成一個小實驗，我們現在有一個數據集，資料集的生成演算法如下：

# -*- coding: utf-8 -*-

from random import random,randint,choice

def hiddenfunction(x,y):
    return x**2 + 2*y + 3*x + 5


def buildhiddenset():
    rows = []
    for i in range(200):
        x=randint(0, 40)
        y=randint(0, 40)
        rows.append([x, y, hiddenfunction(x, y)])
    return rows

if __name__ == '__main__':
    print buildhiddenset()
    

部分資料例如下圖：

視覺化如下，

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from random import random,randint,choice
from mpl_toolkits.mplot3d import axes3d


def hiddenfunction(x,y):
    return x**2 + 2*y + 3*x + 5


def buildhiddenset():
    X = []
    y = []
    for i in range(200):
        x_ = randint(0, 40)
        y_ = randint(0, 40)
        X.append([x_, y_])
        y.append(hiddenfunction(x_, y_))
    return np.array(X), np.array(y)


if __name__ == '__main__':
    # generate a dataset
    X, y = buildhiddenset()

    print "X:", X
    print "y:", y

    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.set_title("3D_Curve")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    # draw the figure, the color is r = read
    figure = ax.plot(X[:, 0], X[:, 1], y, c='r')
    plt.show()

很顯然，必定存在一些函式，可以將（X，Y）對映為結果欄對應的數字，現在問題是這個（些）函式到底是什麼？

從資料分析和數理統計分析的角度來看，這個問題似乎也不是很複雜。我們可以用多元線性迴歸來嘗試擬合數據集。

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from random import random,randint,choice
from mpl_toolkits.mplot3d import axes3d
from sklearn.metrics import accuracy_score
from sklearn.metrics import classification_report, confusion_matrix
import numpy as np
np.set_printoptions(threshold=np.inf)


def hiddenfunction(x,y):
    return x**2 + 2*y + 3*x + 5


def buildhiddenset():
    X = []
    y = []
    for i in range(100):
        x_ = randint(0, 40)
        y_ = randint(0, 40)
        X.append([x_, y_])
        y.append(hiddenfunction(x_, y_))
    return np.array(X), np.array(y)


if __name__ == '__main__':
    # generate a dataset
    X, y = buildhiddenset()

    print "X:", X
    print "y:", y

    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.set_title("3D_Curve")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    # draw the figure, the color is r = read
    #figure = ax.plot(X[:, 0], X[:, 1], y, c='r')
    #plt.show()

    # use Scikit-Learn PolynomialFeature class to constructing parameter terms.
    # a,b,degree=2: [a, b, a^2, ab, b^2]
    # a,b,degree=3: [a, b, a^2, ab, b^2, a^3, a^2b, ab^2, b^3]
    # a,b,c,degree=3: [a, b, c, a^2, ab, ac, b^2, bc, c^2, a^3, a^2b, a^2c, ab^2, ac^2, abc, b^3, b^2c, bc^2, c^3]
    poly_features = PolynomialFeatures(degree=2, include_bias=False)
    # fit the dataset with Polynomial Regression Function, and X_poly is the fitting X result
    X_poly = poly_features.fit_transform(X)
    lin_reg = LinearRegression()
    lin_reg.fit(X_poly, y)
    print(lin_reg.intercept_, lin_reg.coef_)
    y_predict = lin_reg.predict(X_poly)
    y_predict = [int(i) for i in y_predict]
    print "y_predict: ", y_predict
    print "y: ", y
    print 'accuracy for LinearRegression is: {0}'.format(accuracy_score(y, y_predict))
    print 'error for LinearRegression is: {0}'.format(confusion_matrix(y, y_predict))

    # draw the prediction curve
    X_new_1 = np.linspace(0, 40, 100).reshape(100, 1)
    X_new_2 = np.linspace(0, 40, 100).reshape(100, 1)
    X_new = np.hstack((X_new_1, X_new_2))
    # fit the X_new dataset with Polynomial Regression Function, and X_new_poly is the fitting X result
    X_new_poly = poly_features.transform(X_new)
    y_new = lin_reg.predict(X_new_poly)
    y_new = [int(i) for i in y_new]
    #print "X_new: ", X_new
    #print "y_new: ", y_new
    #print "y: ", y

    # draw the prediction line
    figure = ax.plot(X[:, 0], X[:, 1], y, c='r')
    figure = ax.plot(X_new[:, 0], X_new[:, 1], y_new, c='b', linewidth=2, label="Predictions")
    plt.show()

紅色是資料集，藍色線是多項式擬合結果

多項式迴歸的擬合結果如下：

• mean_squared_error for LinearRegression is: 0.55
• accuracy for LinearRegression is: 0.45

因為這是一個迴歸分析問題，因此精度的要求非常高，雖然實際的誤差並不是很高（0.55），使用多項式擬合獲得了0.45的擬合精確度，不算特別高，出錯的點還挺多的。

這麼一看，好像多項式迴歸的效果很差了？並不是這樣的。

我們來對比下模型學習到的多項式引數翻譯成多項式函式，和原始資料集背後的目標函式形式的差距：

print(lin_reg.intercept_, lin_reg.coef_)
(5.0000000000009095, array([ 3.00000000e+00, 2.00000000e+00, 1.00000000e+00, ,
]))
a,b,degree=2: [a, b, a^2, ab, b^2]

模型學習到的：H = 3*X + 2*Y + X^2 + 4.70974246e-17*X*Y -7.81622966e-16*Y^2 + 5    
原始資料集目標函式：H = 3*X + 2*Y + X^2 +  5

可以看到，多項式迴歸出現了一些過擬合的現象，多出了兩項：X*Y、和Y^2，如果沒有這兩項，那麼多項式迴歸得到的方程就是完美的原始方程。

從上面的實驗中，我們可以得到幾點啟發：

• 過擬合問題普遍存在：即使在一個簡單的資料集上，用這麼簡單的多項式迴歸引數，都還是出現了過擬合問題。可見過擬合問題在實際機器學習工程專案中，的確是可能大量存在的。
• 冗餘結構普遍存在：在上面的例子中，X*Y、和Y^2這兩項是多出來的項，模型給這兩個項分配的引數都是一個非常小的數字，小到幾乎可以忽略（分別是4.70974246e-17和7.81622966e-16）。基本上來說，這兩項對最終模型的預測結果的影響是很微小的，如果我們的場景是一個分類任務，那麼這麼小的權重幾乎不會影響到最終的分類預測結果。因此，我們將這種被模型分配權重極小的項，稱為冗餘結構，冗餘結構對分類任務的影響往往非常小。某種程度上來說，冗餘結構緩解了過擬合問題

從這個問題，繼續延伸思考，相信讀者朋友在研究和工程專案中常常會遇到一個有趣的現象：針對某個固定的資料集，我們設計出一個非常精巧的神經網路，擁有上千萬的引數。和另一個人用一個只有十幾萬引數的簡單經典神經網路進行benchmark對比，發現效能並沒有太大的提升，甚至沒有提升。

對這個問題的研究和解釋，已經開始有非常多的學者和科研機構在進行，目前提出了很多的解釋框架（例如teacher-student模型），筆者自己在專案中也同樣遇到了這個問題。一個初步的看法是，對於一個特定場景的資料集來說，是存在一個確定的複雜度上限的。現在如果有兩個模型，一個模型剛剛好達到這個複雜度上限，另一個模型遠遠超過了這個複雜度上限，那麼後面那個模型就會出現所謂的過擬合現象，但是，模型會通過權重分配，將高出複雜度上限之外的多餘項（神經元）都分配0或者非常低的權重，讓這部分神經元稱為冗餘項。最終，不管我們用多麼複雜的模型，起作用的永遠只有那一小部分神經元在發揮作用，預測的效果也一樣好，過擬合問題並沒有影響最終的分類預測結果。

筆者小節：

使用多項式迴歸對一個數據集進行擬合這種做法，本質上是一種先天經驗主義思維，它的科學假設是”複雜可約性“。所謂複雜可約性是複雜理論中的一個概念，指的是一個複雜的系統可以通過一個簡化的模型進行約簡概括，通過歷史的樣本學習這個約簡系統的結構引數，同時假設這個約簡系統能夠比真實的複雜系統推算的更快，從而藉助約簡系統對複雜系統進行未來預測。

上面的話說的有一些繞，簡單來說就是，不管資料集背後的真實目標函式是什麼，我們都可以用像低階多項式函式這種簡單模型來進行擬合學習，並利用學習到的模型對未來可能出現的新資料集進行預測。

這種基於先驗的簡化建模思想在今天的機器學習學術界和工業界應用非常廣泛，也發揮了很好的效果。但其實我們還有另一種看世界的角度，那就是隨機過程。隨機過程事先不做任何假設，而是基於隨機或者遵循某種策略下的隨機，不斷進行自我迭代與優化，從一種混沌狀態逐漸收斂到目標狀態。

0x2：用遺傳程式設計擬合同一個數據集

在這個章節，我們要討論遺傳程式設計，還是上面那個問題，我們現在換一種思維來想：我們所見的資料集背後無非是一個函式公式來決定的，而函式公式又是由一些基本的”數學符號“以及”數學運算子“號組合而成的，數學符號和數學運算子的組合方式我們將其視為一個符號搜尋空間，我們直接去搜索這個符號搜尋空間即可，通過資料集作為反饋，直到收斂為止，即找到了完美的目標函式。

接下來我們逐個小節來分解任務，一步步達到我們的目標。

1. 構造樹狀程式（tree programs）基礎資料結構及操作函式

我們可以用python構造出樹狀程式的基礎資料結構。這棵樹由若干節點組成，根據與之關聯的函式的不同，這些節點又可以擁有一定數量的子節點。

有些節點將會返回傳遞給程式的引數；另一些則會返回常量；還有一些則會返回應用於其子節點之上的操作。

1）fwrapper

一個封裝類，對應於”函式型“節點上的函式。其成員變數包括了函式名稱、函式本身、以及該函式接受的引數個數（子節點）。

class fwrapper:
  def __init__(self,function,params,name):
    self.function=function
    self.childcount=param
    self.name=name

2）node

對應於函式型節點（帶子節點的節點）。我們以一個fwrapper類對其進行初始化。當evaluate被呼叫時，我們會對各個子節點進行求值運算，然後再將函式本身應用於求得的結果。

class node:
  def __init__(self,fw,children):
    self.function=fw.function
    self.name=fw.name
    self.children=children

  def evaluate(self,inp):    
    results=[n.evaluate(inp) for n in self.children]
    return self.function(results)

  def display(self,indent=0):
    print (' '*indent)+self.name
    for c in self.children:
      c.display(indent+1)

3）paramnode

這個類對應的節點只返回傳遞給程式的某個引數。其evaluate方法返回的是由idx指定的引數。

class paramnode:
  def __init__(self,idx):
    self.idx=idx

  def evaluate(self,inp):
    return inp[self.idx]
  
  def display(self,indent=0):
    print '%sp%d' % (' '*indent,self.idx)

4）constnode

返回常量值的節點。其evaluate方法僅返回該類被初始化時所傳入的值。

class constnode:
  def __init__(self,v):
    self.v=v
      
  def evaluate(self,inp):
    return self.v
  
  def display(self,indent=0):
    print '%s%d' % (' '*indent,self.v)

5）節點操作函式

除了基礎數學符號資料結構之外，我們還需要定義一些針對節點的操作函式。

一些簡單的符號運算子（例如add、subtract），可以用lambda內聯方式定義，另外一些稍微複雜的運算子則需要在單獨的語句塊中定義，不論哪種情況，都被會封裝在一個fwrapper類中。

addw=fwrapper(lambda l:l[0]+l[1],2,'add')
subw=fwrapper(lambda l:l[0]-l[1],2,'subtract') 
mulw=fwrapper(lambda l:l[0]*l[1],2,'multiply')

def iffunc(l):
  if l[0]>0: return l[1]
  else: return l[2]
ifw=fwrapper(iffunc,3,'if')

def isgreater(l):
  if l[0]>l[1]: return 1
  else: return 0
gtw=fwrapper(isgreater,2,'isgreater')

flist=[addw,mulw,ifw,gtw,subw]

現在，我們可以利用前面建立的節點類來構造一個程式樹了，我們來嘗試寫一個符號推理程式，

# -*- coding: utf-8 -*-

import gp

def exampletree():
    # if arg[0] > 3:
    #   return arg[1] + 5
    # else:
    #   return arg[1] - 2
    return gp.node(
        gp.ifw, [
            gp.node(gp.gtw, [gp.paramnode(0), gp.constnode(3)]),
            gp.node(gp.addw, [gp.paramnode(1), gp.constnode(5)]),
            gp.node(gp.subw, [gp.paramnode(1), gp.constnode(2)])
        ]
    )

if __name__ == '__main__':
    exampletree = exampletree()

    # expected result = 1
    print exampletree.evaluate([2, 3])

    # expected result = 8
    print exampletree.evaluate([5, 3])

至此，我們已經成功在python中構造出了一個以樹為基礎的語言和直譯器。 

2. 初始化一個隨機種群（函式初始化）

現在我們已經有能力進行形式化符號程式設計了，回到我們的目標，生成一個能夠擬合數據集的函式。首先第一步是需要隨機初始化一個符號函式，即初始化種群。

建立一個隨機程式的步驟包括：

• 建立根節點併為其隨機指定一個關聯函式，然後再隨機建立儘可能多的子節點
• 遞迴地，父節點建立的子節點也可能會有它們自己的隨機關聯子節點

def makerandomtree(pc,maxdepth=4,fpr=0.5,ppr=0.6):
  if random()<fpr and maxdepth>0:
    f=choice(flist)
    children=[makerandomtree(pc,maxdepth-1,fpr,ppr) 
              for i in range(f.childcount)]
    return node(f,children)
  elif random()<ppr:
    return paramnode(randint(0,pc-1))
  else:
    return constnode(randint(0,10))

該函式首先建立了一個節點併為其隨機選了一個函式，然後它遍歷了隨機選中的函式所需的子節點，針對每一個子節點，函式通過遞迴呼叫makerandomtree來建立新的節點。通過這樣的方式，一顆完整的樹就被構造出來了。

僅當被隨機選中的函式不再要求新的子節點時（即如果函式返回的是一個常量或輸入引數時），向下建立分支的過程才會結束。

3. 衡量種群個體的好壞

按照遺傳程式設計演算法的要求，每一輪迭代中都要對種群個體進行定量評估，得到一個個體適應性的排序。

與優化技術一樣，我摩恩必須找到一種衡量題解優劣程度的方法，很多場景下，優劣程度並不容易定量評估（例如網路安全中常常是非黑即白的二分類）。但是在本例中，我們是在一個數值型結果的基礎上對程式進行測試，因此可以很容易通過絕對值誤差進行評估。

def scorefunction(tree,s):
  dif=0
  for data in s:
    v=tree.evaluate([data[0],data[1]])
    dif+=abs(v-data[2])
  return dif

我們來試一試初始化的隨機種群的適應性評估結果，

# -*- coding: utf-8 -*-

import gp

if __name__ == '__main__':
    hiddenset = gp.buildhiddenset()

    random1 = gp.makerandomtree(2)
    random2 = gp.makerandomtree(2)

    print gp.scorefunction(random1, hiddenset)
    print gp.scorefunction(random2, hiddenset)

隨機初始化的函式種群的適應性並不是很好，這符合我們的預期。 

4. 對程式進行變異

當表現最好的程式被選定之後，它們就會被複制並修改以進入到下一代。前面說到，遺傳變異有兩種方式，mutation和crossover，  

1）mutation

變異的做法是對某個程式進行少量的修改，一個樹狀程式可以有多種修改方式，包括：

• 改變節點上的函式
• 改變節點的分支
  • 改變節點所需子節點數目
  • 刪除舊分支
  • 增加新的分支
  • 用全新的樹來替換某一子樹

需要注意的是，變異的次數不宜過多（基因突變不能太頻繁）。例如，我們不宜對整棵樹上的大多數節點都實施變異，相反，我們可以位任何需要進行修改的節點定義一個相對較小的概率。從樹的根節點開始，如果每次生成的隨機數小於該概率值，就以如上所述的某種方式對節點進行變異。

def mutate(t,pc,probchange=0.1):
  if random()<probchange:
    return makerandomtree(pc)
  else:
    result=deepcopy(t)
    if hasattr(t,"children"):
      result.children=[mutate(c,pc,probchange) for c in t.children]
    return result

2）crossover

除了變異，另一種修改程式的方法被稱為交叉或配對，即：從本輪種群中的優秀適應著中，選出兩個將其進行部分子樹交換。執行交叉操作的函式以兩棵樹作為輸入，並同時開始向下遍歷，當到達某個隨機選定的閾值時，該函式便會返回前一棵樹的一份拷貝，樹上的某個分支會被後一棵樹上的一個分支所取代。通過同時對兩棵樹的即時遍歷，函式會在每棵樹上大致位於相同層次的節點處實施交叉操作。

def crossover(t1,t2,probswap=0.7,top=1):
  if random()<probswap and not top:
    return deepcopy(t2) 
  else:
    result=deepcopy(t1)
    if hasattr(t1,'children') and hasattr(t2,'children'):
      result.children=[crossover(c,choice(t2.children),probswap,0) 
                       for c in t1.children]
    return result

讀者朋友可能會注意到，對於某次具體的變異或者交叉來說，新的種群個體並不一定會帶來更好的效能，實際上，新種群個體的效能幾乎完全是隨機的。從生物進化論的角度來說，遺傳變異是無方向的，隨機的，遺傳變異的目標僅僅是引入多樣性，造成演化的是環境選擇壓（資料集的誤差反饋）。

5. 持續迭代演化

現在，我們將上面的步驟串起來，讓遺傳演化不斷的迴圈進行。本質上，我們的思路是要生成一組隨機程式並擇優複製和修改，然後一直重複這一過程直到終止條件滿足為止。

def getrankfunction(dataset):
  def rankfunction(population):
    scores=[(scorefunction(t,dataset),t) for t in population]
    scores.sort()
    return scores
  return rankfunction
  
    

def evolve(pc,popsize,rankfunction,maxgen=500,
           mutationrate=0.1,breedingrate=0.4,pexp=0.7,pnew=0.05):
  # Returns a random number, tending towards lower numbers. The lower pexp
  # is, more lower numbers you will get
  def selectindex():
    return int(log(random())/log(pexp))

  # Create a random initial population
  population=[makerandomtree(pc) for i in range(popsize)]
  for i in range(maxgen):
    scores=rankfunction(population)
    print "function score: ", scores[0][0]
    if scores[0][0]==0: break
    
    # The two best always make it
    newpop=[scores[0][1],scores[1][1]]
    
    # Build the next generation
    while len(newpop)<popsize:
      if random()>pnew:
        newpop.append(mutate(
                      crossover(scores[selectindex()][1],
                                 scores[selectindex()][1],
                                probswap=breedingrate),
                        pc,probchange=mutationrate))
      else:
        # Add a random node to mix things up
        newpop.append(makerandomtree(pc))
        
    population=newpop
  scores[0][1].display()    
  return scores[0][1]

上述函式首先建立一個隨機種群，然後迴圈至多maxgen次，每次迴圈都會呼叫rankfunction對程式按表現從優到劣的順序進行排列。表現優者會不加修改地自動進入到下一代，我們稱這樣的方法為精英選拔髮（elitism）。

至於下一代中的其他程式，則是通過隨機選擇排名靠前者，再經過交叉和變異之後得到的。

這個過程是一直重複下去，知道某個程式達到了完美的擬合適配（損失為0），或者重複次數達到了maxgen次為止。

evolve函式有多個引數，用以從不同方面對競爭環境加以控制，說明如下：

• rankfunction：對應一個函式，將一組程式從優到劣的順序進行排列
• mutationrate：代表發生變異的概率
• breedingrate：代表發生交叉的概率
• popsize：初始種群的大小
• probexp：表示在構造新種群時，”選擇評價較低的程式“這一概率的遞減比例。該值越大，相應的篩選過程就越嚴格，即只選擇評價最高的多少比例的個體作為複製物件
• probnew：表示在構造新種群時，”引入一個全新的隨機程式“的概率，該引數和probexp是”種群多樣性“的重要決定引數

# -*- coding: utf-8 -*-

import gp

if __name__ == '__main__':
    rf = gp.getrankfunction(gp.buildhiddenset())
    gp.evolve(2, 500, rf, mutationrate=0.2, breedingrate=0.1)

程式執行的非常慢，在筆者的mac上運行了15min才最終收斂到0。有意思的是，儘管這裡給出的解是完全正確的，但是它明顯比我們資料集背後的真實目標函式要複雜得多，也就是說發生了過擬合。

但是，我們如果運用一些代數知識，將遺傳程式設計得到函式進行約簡，會發現它和目標函式其實是等價的（p0為X，p1為Y）。

((X+6)+Y)+X + (if( (X*Y)>0 ){X}else{X} + X*X) + (X + (Y - (X + if(6>0){1}else{0})) )
# X*Y恆大於0
2*X + Y + 6 + X + X**2 + (X + (Y - (X + if(6>0){1}else{0})) )
# 6恆大於0
2*X + Y + 6 + X + X**2 + X + Y - X + 1
X**2 + 3*X + 2*Y + 5

可以看到，遺傳程式設計這種基於內顯性的構造方式，可以在形式上得到一個全域性最優解，這點上比基於優化演算法的逼近方法要好。

同時，上述例子告訴我們遺傳程式設計的一個重要特徵：遺傳演算法找到的題解也許是完全正確的，亦或是非常不錯的。但是通常這些題解遠比真實的目標函式要複雜得多。在遺傳程式設計得到的題解中，我們發現有很多部分是不做任何工作的，或者對應的是形式複雜，但始終都只返回同一結果的公式，例如"if(6>0){1}else{0})"，只是1的一種多餘的表達方式而已。

0x3：從遺傳程式設計優化結果看過擬合問題的本質

從上一章的遺傳優化結果我們可以看到，在資料集是充分典型集的情況下，過擬合是不影響模型的收斂的。

遺傳程式設計的這種冗餘性和優化演算法和深度神經網路中的冗餘結構本質上是一致的，這是一個冗餘的過擬合現象，即程式的題解非常複雜，但是卻對最終的決策沒有影響，唯一的缺點就是浪費了很多cpu時鐘。

但是另一方面，這種冗餘過擬合帶來的一個額外的風險是，”如果資料集是非典型的，那麼過擬合就會導致嚴重的後果“。

我們需要明白的是，過擬合的罪魁禍首不是模型和優化演算法，而恰恰是資料集本身。在本例中我們清楚地看到，當我們能夠得到一個完整的典型集訓練資料時，過擬合問題就會退化成一個冗餘魯棒可約結構。

但是反之，如果我們的資料集因為系統噪聲或採樣不完全等原因，沒有拿到目標函式的典型集，那麼由於複雜模型帶來的過擬合問題就會引發很嚴重的預測偏差。我們來稍微修改一下程式碼，將原始資料集中隨機剔除1/10的資料，使資料的充分性和典型性下降，然後再來看遺傳程式設計最後的函式優化結果，

# -*- coding: utf-8 -*-

import gp

if __name__ == '__main__':
    hiddenset = gp.buildhiddenset()
    # 按照 5% 模來取樣，即剔除1/10的資料，模擬取樣不完全的情況
    cnt = 0
    for i in hiddenset:
        if cnt % 10 == 0:
            hiddenset.remove(i)
        cnt += 1
    print hiddenset

    rf = gp.getrankfunction(hiddenset)
    gp.evolve(2, 500, rf, mutationrate=0.2, breedingrate=0.1)

得到的函式約簡結果為： 

if( ((Y+4) + if(Y>X){1}else{0} ) > 0){1}else{0} + (Y+4) + X^2 + (Y + (X + (X+X)))
# if( ((Y+4) + if(Y>X){1}else{0} ) > 0){1}else{0} 不可約
if( ((Y+4) + if(Y>X){1}else{0} ) > 0){1}else{0} + X**2 + 3*X + 2*Y + 4

# 真實的目標函式為:
X**2 + 3*X + 2*Y + 5

可以看到，在資料集不完整的情況下，遺傳演算法就算完美擬合了訓練集，但是也無法真正逼近目標函式，但這不是遺傳演算法的問題，而是受到資料集的制約。

更重要的是，因為資料集的非典型性（資料概率分佈缺失），導致模型引入了真正的”過擬合複雜結構“，即”if( ((Y+4) + if(Y>X){1}else{0} ) > 0){1}else{0}“，這是一個區間函式。要知道，這僅僅是一個非常小的例子，尚且引入瞭如此的不確定性，在更復雜和更復雜的問題中，資料集的概率分佈缺失會引發更大的”多餘過擬合複雜結構問題“，影響程度的多少，根據資料缺失的程度而定。

這反過來提醒了我們這些資料科學工作者，在解決一個實際問題的時候，不要過分糾結你的模型是不是足夠好，層數是不是足夠深，而要更多地關注你的資料集，資料集的質量直接決定了最終的模型效果。更進一步地說，如果你能找到一種方法，能100%拿到目標函式的資料全集，那麼恭喜你，隨便用一個機器學習模型都可以取得很好的效果。

0x4：多樣性的重要性

對於遺傳程式設計，我們還需要再談一個關鍵問題，即多樣性問題。

我們看到evolve函式中，會將最優的個體直接進入下一代，除此之外，對排名之後的個體也會按照比例和概率選擇性地進行復制和修改以形成新的種群，這種做法有什麼意義呢？

最大的問題在於，僅僅選擇表現最優異的少數個體，很快就會使種群變得極端同質化（homogeneous），或稱為近親交配。儘管種群中所包含的題解，表現都非常不錯，但是它們彼此間不會有太大的差異，因為在這些題解間進行的交叉操作最終會導致群內的題解變得越來越相似。我們稱這一現象為達到區域性最大化（local maxima）。

對於種群而言，區域性最大化是一種不錯的狀態（即收斂了），但還稱不上最佳的狀態。因為處於這種狀態的種群裡，任何細小的變化都不會對最終的結果產生太大的變化。這就是一個哲學上的矛盾與對立，收斂穩定與發散變化是彼此對立又統一的，完全偏向任何一方都是不對的。

事實表明，將表現極為優異的題解和大量成績尚可的題解組合在一起，往往能夠得到更好的結果。基於這個原因，evolve提供了兩個額外的引數，允許我們對篩選程序中的多樣性進行調整。

• 通過降低probexp的值，我們允許表現較差的題解進入最終的種群之中，從而將”適者生存（survival of fittest）“的篩選過程調整為”最適應者及其最幸運者生存（survival of the fittest and luckiest）“
• 通過增加probnew的值，我們還允許全新的程式被隨機地加入到種群中

這兩個引數都會有效地增加進化過程中的多樣性，同時又不會對程序有過多的擾亂，因為，表現最差的程式最終總是會被剔除掉的（遺傳程式設計的馬爾科夫收斂性）。

Relevant Link:  

《集體智慧程式設計》 

5. 用遺傳程式設計自動得到正則表示式生成器 - Regex Golf Problem

0x1：問題描述

我們需要生成一段正則表示式，這個正則表示式需要能夠匹配到所有的M資料集，同時不匹配所有的U資料集，且同時還要儘量短，即不能是簡單的M資料集的並集拼接。

定義一個目標（損失）函式來評估每次題解的好壞，

，其中n_M代表匹配M的個數，n_U代表匹配U的個數，w_I代表獎勵權重，r代表該正則表示式的長度

演算法優化的目標是使上式儘量大。 

0x2：從一個貪婪演算法說起 

在討論遺傳程式設計之前，我們先從常規思路，用一種貪婪迭代演算法來解決這個問題。我們的資料集如下，

# -*- coding: utf-8 -*-

from __future__ import division
import re
import itertools

def words(text):
    return set(text.split())

if __name__ == '__main__':
    M = words('''afoot catfoot dogfoot fanfoot foody foolery foolish fooster footage
        foothot footle footpad footway hotfoot jawfoot mafoo nonfood padfoot prefool sfoot unfool''')

    U = words('''Atlas Aymoro Iberic Mahran Ormazd Silipan altared chandoo crenel crooked
        fardo folksy forest hebamic idgah manlike marly palazzi sixfold tarrock unfold''')

    print M & U

首先，確認了U和M之間不存在交集，這道題理論上是有解的，否則無解。

1. 定義可行解判斷條件

我們先準確定義出什麼時候意味著得到了一個可行解，

# -*- coding: utf-8 -*-

from __future__ import division
import re
import itertools


def words(text):
    return set(text.split())


def mistakes(regex, M, U):
    "The set of mistakes made by this regex in classifying M and U."
    return ({"Should have matched: " + W for W in M if not re.search(regex, W)} |
            {"Should not have matched: " + L for L in U if re.search(regex, L)})


def verify(regex, M, U):
    assert not mistakes(regex, M, U)
    return True


if __name__ == '__main__':
    M = words('''afoot catfoot dogfoot fanfoot foody foolery foolish fooster footage
        foothot footle footpad footway hotfoot jawfoot mafoo nonfood padfoot prefool sfoot unfool''')

    U = words('''Atlas Aymoro Iberic Mahran Ormazd Silipan altared chandoo crenel crooked
        fardo folksy forest hebamic idgah manlike marly palazzi sixfold tarrock unfold''')

    some_answer = "a*"

    print mistakes(some_answer, M, U)

可以看到，當我們輸入正則”a*“的時候出現了很多錯誤，顯然”a*“不是我們要的答案。讀者朋友可以試著輸入”foo“試試。

2. 尋找可行解的策略

• 對M中的每個詞（短語）都進行一次正則候選集構造，包括以下步驟：
  • 遍歷M中的每一個詞的每一次字元，並過濾掉特殊字元（*+?^$.[](){}|\\），然後在中間遍歷插入”*+？“，這是字元級別的混合交叉
  • 對M中的每一個詞都加上首尾定界符，例如”^it$“，得到一個wholes詞集
  • 對wholes詞集進行ngram切分，得到一個ngram詞集，例如對於詞”^it$“來說，可以得到{'^', 'i', 't', '$', '^i', 'it', 't$', '^it', 'it$', '^it$'}，作為一個正則串池。可以這麼理解，這個池中的每個正則串都至少能命中一個M中的元素
  • 遍歷上一步得到的正則串池中所有元素，逐字元用”.“字元進行替換，例如對於"^it$"來說，可以得到{'^it$', '^i.$', '^.t$', '^..$'}
  • 遍歷上一步dotify的詞集，逐字元遍歷插入”*+？“這種repetition控制符，例如對於”a.c“來說，可以得到{'a+.c', 'a*.c', 'a?.c','a.c+', 'a.*c', 'a.?c','a.+c', 'a.c*', 'a.c?'}，需要注意的是，在首位定界符前後不要加repetition控制符，同時不要同時加入2個repetition控制符
• 從題解候選集中，篩選出至少能夠匹配一個以上M，但是不匹配U的正則子串，這一步得到一個題解正則候選子集。這是一個貪婪迭代式的思想，它不求一步得到一條能夠匹配所有M的正則，而是尋找一些能夠解決一部分問題的正則子串，將困難問題分而治之
• 使用OR拼接將題解正則候選子集拼接起來，例如”ab | cd“

上面構造正則候選集的過程說的可能有些抽象，這裡通過程式碼示例來說明一下，

# -*- coding: utf-8 -*-

from __future__ import division
import re
import itertools


OR  = '|'.join # Join a sequence of strings with '|' between them
cat = ''.join  # Join a sequence of strings with nothing between them
Set = frozenset # Data will be frozensets, so they can't be mutated.


def words(text):
    return set(text.split())


def mistakes(regex, M, U):
    "The set of mistakes made by this regex in classifying M and U."
    return ({"Should have matched: " + W for W in M if not re.search(regex, W)} |
            {"Should not have matched: " + L for L in U if re.search(regex, L)})


def verify(regex, M, U):
    assert not mistakes(regex, M, U)
    return True


def matches(regex, strings):
    "Return a set of all the strings that are matched by regex."
    return {s for s in strings if re.search(regex, s)}


def regex_parts(M, U):
    "Return parts that match at least one winner, but no loser."
    wholes = {'^' + w + '$' for w in M}
    parts = {d for w in wholes for p in subparts(w) for d in dotify(p)}
    return wholes | {p for p in parts if not matches(p, U)}


def subparts(word, N=4):
    "Return a set of