論文式程式設計

meetrice發表於2016-02-06

文學程式設計

文學程式設計(Literate programming)的一些概念,上個世紀 70 年代就有人提出來了。

文學程式設計的思想非常簡單,就是將那些為了能被編譯器/直譯器正確識別而編寫的程式碼打碎,然後用人類語言將它們編織到文件中,這種文件就是文學程式設計的原始檔。這一概念第一次被完整的實現,是 Knuth 開發的 WEB 工具(此 WEB 並非現代漫天飛舞的那個 Web)。Knuth 的神作——TeX 系統便是藉助 WEB 開發的。

WEB 工具由 tangle 與 weave 這兩個程式構成。tangle 程式從文學程式設計的原始檔中提取複合編譯器/直譯器邏輯的程式程式碼。weave 程式將文學程式設計的原始檔轉換為 TeX 原始檔,然後由 TeX 系統排版處理,生成程式文件。這種程式文件使作者能在以後的任何時間重新找到自己的思路,也能使其他程式設計師更容易理解程式的建構過程。

在程式碼裡寫大量的註釋,或者像 Doxygen 之類從具有特定註釋格式的程式碼中產生文件的工具,這些都不算文學程式設計。文學程式設計強調的是,程式碼出現的順序應該按照人的邏輯,而不是編譯器的。

很近的遙遠

文學程式設計,從未真正的被廣泛應用,但是每個程式猿都曾經看到過它的影子。閱讀一本講 Linux 核心的書籍,遠比閱讀 Linux 核心原始碼容易的多。

每一本講程式設計技術的書籍都是以類似『文學程式設計』的方式寫成的。我們在閱讀這些書籍時,總是先關心作者說了些什麼,然後再看他給出的示例程式碼。事實上,有本很古老但是很有名而且現在也不過時的書,它的中文譯本叫《C 語言介面與實現——建立可重用軟體的技術》,這本書就是基於文學程式設計的方式實現的,作者用的文學程式設計工具叫 noweb——後面我會介紹它。

既然文學程式設計如此的美好,那麼為何它一直沒有被廣泛應用呢?答案很簡單,寫書(不是那種爛書)要比寫程式難得多。因此,寫面向人類讀者的文學程式要比寫面向機器的程式碼難得多。noweb 的作者在 noweb 的原始檔裡說:One of my observations is that the cost of creating a high quality, well-documented literate program adds 1-3 times of amount of effort it took to create the program in the first place.

大部分程式猿並不那麼熱愛程式設計,他們僅僅是將程式設計當作養家餬口的一種手段。採用文學程式設計方式幹活,要比面向機器程式設計多付出 1 到 3 倍的努力,而薪酬卻不變……太多的正常人是絕對不會這麼虐待自己的。

這個世界從不缺乏足夠好的東西,只是缺乏足夠好的人。

它文學麼?

什麼是文學?據說,文學是以語言文字為工具,形象化地反映客觀現實、表現作家心靈世界的藝術。

雖然文學程式設計用的都是語言文字,但是它的目標卻不是形象化的反映客觀現實或表現作者心靈世界的,它的目標是寫出讓人類與機器都能能容易的解讀的『文件』。程式設計的本質,是讓計算機執行確定的計算。計算的本質是數學意義上的。如果數學還沒有變成文學,那麼程式設計永遠也無法變成文學。

這一點決定了文學程式設計並不能用來搞文學,所以不要被它的名字所迷惑。要搞文學創作,只能熱愛生活,認真思考,錘鍊文字。

論文式程式設計

Knuth 明確提出了文學程式設計的概念,並付諸於實踐,開發了 TeX、MetaFont 以及 MMIX 元模擬器這些大的程式,此外還出版了一本闡述文學程式設計的專著。文學程式設計對於他的重要作用,用他自己的話來說,就是:『文學程式設計確實是由TeX專案衍生出來的最重要的東西。它不僅讓我前所未有地更快地寫和維護可靠性更高的程式,而且成為我自20世紀80年代以來的最大的快樂之源——它有時實際上是不可或缺的。我做的其它一些大程式,比如 MMIX 元模擬器,用我見過的任何一種其它的方法論是無法寫出來的。其複雜性讓我有限的智慧望而卻步。沒有文學程式設計,我的整個事業規劃就會轟然倒塌。……文學程式設計是你更上一層樓的必要工具。』

如果我們打算用文學程式設計的方式程式設計,那麼必須要注意,Knuth 是計算機領域的科學家!我不是在鼓吹他的個人頭銜,而是強調他的身份。科學家,是這個世界上最會寫論文的一群人。他們探索未知區域,在失敗中前進,認真總結自己的發現,最後以嚴肅的論文或專著的形式將自己的發現公佈於世。TeX 對於我們來說是一個程式,但它對於 Knuth 來說是一本講述計算機排版技術的專著。也就是說,如果想用好文學程式設計,那麼首先你得學會如何寫文章以及如何寫科技文獻。

寫一個程式與寫一本專著,寫一個程式模組與寫一篇科技論文,它們之間是存在著對應關係的。寫程式,首先要明確程式所解決的問題,然後將問題拆分為更小的問題,每個問題都用一個比較小的程式模組來解決,最後將各個模組組裝起來得到程式。做科研的人,首先要明確自己要研究的主題,然後將主題分解為一些小的主題,對小的主題展開研究,每個研究成果都以論文的形式發表,最後將整個主題的相關論文整合起來形成專著。

所以,我覺得文學程式設計,應該叫論文式程式設計,至少也該叫『文式程式設計』。如果按後者來理解,那麼本文的標題應該這麼讀『論(文式程式設計)』:)

如果我們像寫一篇論文那樣來寫一個程式模組,大致的過程應該是:

  • 引言(背景):這個程式模組要解決什麼問題?要解決這個問題,有哪些現有的資源?我該如何利用這些資源來解決問題?

  • 方法(演算法與程式):準確描述利用現有資源解決問題的全部步驟,有序的組織程式程式碼。

  • 結果(單元測試):驗證方法是否正確。

  • 討論:使用這個模組需要注意什麼,它還存在哪些不足,應該怎樣彌補。

其實,我們在寫每一個程式模組時,都會經歷上述的思考與實現過程,我們最終所得到的是程式碼,這個過程中我們的那些思考活動卻很少被記錄下來。很多人說,我要去讀 XXXX 專案的原始碼,這是學習 XXXX 的最好方法。這種觀點其實並不正確,因為你閱讀原始碼的過程,實際上就是在猜測這些程式碼當初是怎樣寫出來的。大可以放心,無論你怎麼猜測,你只能得到一個很模糊的結果,因為那些原本很確切的資訊已經永遠的丟失了,甚至連當初寫這些程式碼的人可能也想不起來了,他們留給你的只是一個巨大的迷宮。

雖然有一些程式碼是自明的,但是,顯然這些程式碼也都是非常簡單的。對一個矩陣進行奇異值分解(SVD)的程式碼,無論怎麼寫,它也無法是自明的,除非你去閱讀一篇闡述矩陣奇異值分解演算法的論文。

示例

作為示例,我要用論文式程式設計的方式來寫一個遺傳演算法的程式。這個程式的原始碼如下:

% -*- mode: Noweb; noweb-code-mode: python-mode -*-
\title{Hello!遺傳演算法}
\cprotect\author{garfileo\\ \verb|lyr_m2@live.cn|}
\date{\today}
\maketitle

\tableofcontents
\newpage

\section*{引言}

這篇文章講述如何利用遺傳演算法解決一個二元函式的最大值求解問題。由於我對遺傳演算法的理解處於菜鳥級別,所以本文所講的方法以及所寫的程式不一定正確。之所以寫這篇文章,是因為我已經煩透了教科書或論文裡對遺傳演算法那麼刻板的敘述,所以很想寫一篇稍微輕鬆一點的入門文件,娛樂一下。

\section{問題}

這個二元函式是這樣的:

$$f(x,\,y)=0.5-\frac{\sin^2{\sqrt{x^2+y^2}-0.5}}{1+0.001(x^2+y^2))^2}$$

要是我能夠在大腦中直接生成這個函式的影象就好了,可惜我不能夠,所以用 gnuplot 畫了一下。

\begin{figure}[htbp]
\centering
\includegraphics[width=6cm]{f.png}
\caption[目標函式]{待求最大值的目標函式}
\end{figure}

這個函式像是平靜的池塘裡丟了一顆小石子激起的波紋。我們的任務是計算它在 $x\in [-10,\,10],\;y\in [-10,\,10]$ 範圍之內的最大值。

這個函式有無限個極大值,但是僅有一個最大值,位於 $(0, 0)$ 點,值為 1。如果你的微積分學知識還沒有遺忘,可以用數學方法求解一番。不幸的是,我已經忘光了,所以我只好用遺傳演算法進行求解。遺傳演算法的特點之一就是:{\bf 不需要求導或其他輔助知識,而只需要影響一些可以影響搜尋方向的目標函式和相應的適應度函式}。所謂目標函式,就是要求解的函式,也就是上述的那個函式。至於適應度函式,下文再行介紹。

\section{建立染色體}

我唯一接觸到生物學是在我的初中時代。就讀的那個初中學校是一個落後的鄉村中學,不過卻擁有一個很好的教生物的老師,但是悲劇的是我在那時是一個不喜歡上課的懵懂無知的少年。現在為了理解遺傳演算法,我只好將『染色體』理解成一根帶子,上面寫著一組資料。據說這組資料記錄著我們應該長成什麼樣子,具備什麼樣的天賦,可能會生什麼疾病等內容。如果上帝能夠將『語言程式』記錄在我們的染色體中,也許我們剛生下來就可以說上百種人類語言還有火星語了。

雖然我們不是上帝,但是我們也可以創造染色體,例如 $000110001100$ 或者 $000XXX00XXX0X$. 這是一件很容易的事情,而真正困難的是如何在染色體中記錄資訊。由於用二進位制來表示染色體比較方便程式計算,所以本文選擇了這種最簡單的方式。

現在,嘗試為 $f(x, y)$ 的最大值所對應的 $x$ 和 $y$ 的值構造染色體。也就是說,要在一組二進位制位中儲存 $f(x, y)$ 的定義域中的數值資訊。

顯然,函式 $f(x, y)$ 的定義域所包含的數值是無限多的,但是基於取樣的辦法可以得到有限集。例如,對於 $[-10,\,10]$ 這個區間,我們可以將它平均劃分為 $20\times 10^6$ 個子區間,便得到精度為 8 位,小數位為 6 位的一組數值,個數為 $20\times 10^6 + 1$ 。若用一組二進位制位形式的染色體來表示這個數值集合,那麼我們還要考慮所用二進位制位的長度。由於 $2^{24}<20\times 10^6 + 1< 2^{25}$,因此我們可以將染色體長度確定為 25 位,因為只有如此才可以讓足夠多的染色體表示那麼多的數值,同時又不至於太浪費。雖然長度為 25 的二進位制位所能表示的數值個數要多於 $20\times 10^6 + 1$,但是這並沒有負面作用,相反,它可以更精確的表示區間 $[-10,\, 10]$ 中數值。

現在,我們已經建立了一種 25 位長度的二進位制位型別的染色體,那麼對於任意一個這樣的染色體,我們如何將其復原為 $[-10,\,10]$ 這個區間中的數值呢?很簡單,只需要使用下面的公式:

$$f(c) = -10.0 + c\cdot\frac{10.0 - (-10.0)}{2^{25} - 1}$$

例 如 $0000 0000 0000 0000 0000 0000 0000 0$ 和 $1111 1111 1111 1111 1111 1111 1$ 這兩個二進位制數,將其化為 10 進位制數,代入上式,可得 -10.0 和 10.0。這意味著長度為 25 位的二進位制數總是可以通過上式轉化為 $[-10,\,10]$ 區間中的數。

\section{個體、種群與進化}

染色體表達了某種特徵,這種特徵的載體,可以稱為『個體』。例如,我本人就是一個『個體』,我身上載有 23 對染色體,也許我的相貌、性別、性格等因素主要取決於它們。眾多個體便構成『種群』。

對於本文所要解決的二元函式最大值求解問題,個體可採用上一節所構造的染色體表示,並且數量為 2 個,其含義可理解為函式 $f(x, y$) 定義域內的一個點的座標。許多這樣的個體便構成了一個種群,其含義為一個二維點集,包含於對角定點為 $(-10.0, -10.0)$ 和 $(10.0, 10.0)$ 的正方形區域。

也許有這樣一個種群,它所包含的個體對應的函式值會比其他個體更接近於函式 $f(x, y)$ 的理論最大值,但是它一開始的時候可能並不比其他個體優秀,它之所以優秀是因為它選擇了不斷的進化,每一次的進化都要儘量保留種群中的優秀個體,淘汰掉不理想的個體,並且在優秀個體之間進行染色體交叉,有些個體還可能出現變異。種群的每一次進化後,必定會產生一個最優秀的個體。種群所有世代中的那個最優個體也許就是函式 $f(x, y)$ 的最大值對應的定義域中的點。如果種群不休止的進化,它總是能夠找到最好的解。但是,由於我們的時間是有限的,有可能等不及種群的最優進化結果,通常是在得到了一個看上去還不錯的解時,便終止了種群的進化。

那麼,對於一個給定的種群,如何賦予它進化的能力呢?

\begin{itemize}
\item {\bf 選擇}:對於種群的每一代個體,可以用一個適應度函式(也叫評估函式)計算個體的適應度,根據適應度可以計算出個體的生存機率。適應度較大的個體被保留的可能性也較大,反之被淘汰的可能性較大。
\item {\bf 交叉}:在一定的概率下對兩個個體的染色體進行交叉重組,從而得到兩個新個體。
\item {\bf 變異}:些個體的染色體會以一定的概率發生變化。
\end{itemize}

達爾文的進化論也許並不正確,但是它對於我們運用這種理論來計算問題並沒有什麼錯誤的影響。我們不管人類是否是由猿猴進化來的,還是由別的什麼生物。那些進化論的反對者總是想用自己的理論推翻進化論,不過他們的理論卻往往無法用於計算!基督徒們相信世界末日,也許只是因為上帝的時間也很有限,等不及人類進化到最優解,於是就設定了人類進化的最大世代數。

\section{種群}

如果你不熟悉 python 語言,那麼請原諒我使用了它。我將種群宣告為 python 的一個類:

<<種群>>= 
class Population
@

種群的初始化過程就是 `Population` 類的初始化函式:

<<種群初始化>>=
def __init__ (self, size, chrom_size, cp, mp, gen_max):
    self.individuals = []          # 個體集合
    self.fitness = []              # 個體適應度集合
    self.selector_probability = [] # 個體選擇概率集合
    self.new_individuals = []      # 新一代個體集合
    
    self.elitist = {'chromosome':[0, 0],
                    'fitness':0,
                    'age':0}       # 最佳個體的資訊
    
    self.size = size # 種群所包含的個體數
    self.chromosome_size = chrom_size # 個體的染色體長度
    self.crossover_probability = cp   # 個體之間的交叉概率
    self.mutation_probability = mp    # 個體之間的變異概率
    
    self.generation_max = gen_max # 種群進化的最大世代數
    self.age = 0                  # 種群當前所處世代
     
    # 隨機產生初始個體集,並將新一代個體、適應度、選擇概率等集合以 0 值進行初始化
    v = 2 ** self.chromosome_size - 1
    for i in range (self.size):
        self.individuals.append ([random.randint (0, v), random.randint (0, v)])
        self.new_individuals.append ([0, 0])
        self.fitness.append (0)
        self.selector_probability.append (0)
@

程式碼中的 [[self]] 就是種群的例項,下文中也是如此。

\section{選擇}

可以簡單的模擬出『物競天擇』的效果:將種群的各個個體擺在一個輪盤上,然後轉一下輪盤,將盤外的指標所指向的個體保留下來,然後接著轉輪盤,接著選擇,直至產生一組與種群原有個體數量一致的個體,這就是我們所選擇的下一代。這種賭博不違法。

要模擬這個輪盤賭博機制,首先需要構造個體適應度評價機制:

<<物競天擇機制>>=
def decode (self, interval, chromosome):
    d = interval[1] - interval[0]
    n = float (2 ** self.chromosome_size -1)
    return (interval[0] + chromosome * d / n)

def fitness_func (self, chrom1, chrom2):
    interval = [-10.0, 10.0]
    (x, y) = (self.decode (interval, chrom1), 
              self.decode (interval, chrom2))
    n = lambda x, y: math.sin (math.sqrt (x*x + y*y)) ** 2 - 0.5
    d = lambda x, y: (1 + 0.001 * (x*x + y*y)) ** 2
    func = lambda x, y: 0.5 - n (x, y)/d (x, y)
    return func (x, y)
    
def evaluate (self):
    sp = self.selector_probability
    for i in range (self.size):
        self.fitness[i] = self.fitness_func (self.individuals[i][0], 
                                             self.individuals[i][1])
    ft_sum = sum (self.fitness)
    for i in range (self.size):
        sp[i] = self.fitness[i] / float (ft_sum)
    for i in range (1, self.size):
        sp[i] = sp[i] + sp[i-1]
@

[[decode]] 函式可以將染色體 [[chromosome]] 對映為區間 [[interval]] 之內的數值。[[fitness_func]] 是適應度函式,可以根據個體的兩個染色體計算出該個體的適應度,這裡直接採用了本文所要求解的目標函式

$$f(x,\,y)=0.5-\frac{\sin^2{\sqrt{x^2+y^2}-0.5}}{1+0.001(x^2+y^2))^2}$$

作為適應度函式。

[[evaluate]] 函式用於評估種群中的個體集合 [[self.individuals]] 中各個個體的適應度,即將各個個體的 2 個染色體代入 [[fitness_func]] 函式,並將計算結果儲存在 [[self.fitness]] 列表中,然後將 [[self.fitness]] 中的各個個體適應度除以所有個體適應度之和,得到各個個體的生存概率。為了適合輪盤賭博遊戲,需要將個體的生存概率進行疊加,從而計算出各個個體的選擇概率。例如有 5 個個體,根據其適應度計算的生存概率與選擇概率如表 \ref{table:選擇概率計算示例} 所示。

\begin{table}[H]
\centering
\caption{選擇概率的計算結果示例}
\label{table:選擇概率計算示例}
\begin{tabular}{cccc}
\toprule
\bf 個體 & \bf 適應度 & \bf 生存概率 & \bf 選擇概率 \\
\midrule
1 & 0.9042845033795694 & 0.28693981857759787 & 0.28693981857759787 \\
2 & 0.5588628304907922 & 0.17733356990137467 & 0.46427338847897254 \\
3 & 0.6899948769706024 & 0.21894326849291637 & 0.6832166569718889 \\
4 & 0.3114709778723004 & 0.09883330472749545 & 0.7820499616993843 \\
5 & 0.6868647339474463 & 0.21795003830061557 & 0.9999999999999999 \\
\bottomrule
\end{tabular}
\end{table}

有了這些資料,便可以構造圖 \ref{fig:輪盤} 所示的輪盤賭博機了。

\begin{figure}[h]
\centering
\includegraphics[width=4cm]{selector.png}
\caption{輪盤賭博機}
\label{fig:輪盤}
\end{figure}

這樣的輪盤賭博機,可用 python 程式碼表示為:

<<物競天擇機制>>=
def select (self):
    (t, i) = (random.random (), 0)
    for p in self.selector_probability:
        if p > t:
            break
        i = i + 1
    return i
@ 

\section{染色體交叉模擬}

<<染色體交叉機制>>=
def cross (self, chrom1, chrom2):
    p = random.random ()
    n = 2 ** self.chromosome_size -1
    if chrom1 != chrom2 and p < self.crossover_probability:
        t = random.randint (1, self.chromosome_size - 1)
        mask = n << t
        (r1, r2) = (chrom1 & mask, chrom2 & mask)
        mask = n >> (self.chromosome_size - t)
        (l1, l2) = (chrom1 & mask, chrom2 & mask)
        (chrom1, chrom2) = (r1 + l2, r2 + l1)
    return (chrom1, chrom2)
@

[[cross]] 函式可以將兩個染色體進行交叉配對,從而生成 2 個新染色體。

此處使用染色體交叉方法很簡單,先生成一個隨機概率 [[p]],如果兩個待交叉的染色體不同並且 [[p]] 小於種群個體之間的交叉概率 [[self.crossover_probability]],那麼就在 $[0, \text{self.chromosome\_size}]$ 中間隨機選取一個位置,將兩個染色體分別斷為 2 截,然後彼此交換一下。例如:

\begin{verbatim}
1000 1101 1100 0010 0001 0110 1
0001 0011 1111 1001 0010 1110 0
\end{verbatim}

\noindent 在第 10 位處交叉,結果為:

\begin{verbatim}
1000 1101 1100 0011 0010 1110 0
0001 0011 1111 1000 0001 0110 1
\end{verbatim}

這種染色體交叉方法叫做{\bf 單點交叉}。如果不嫌麻煩,也可以使用{\bf 多點交叉}。

\section{染色體變異}

<<染色體變異機制>>=
def mutate (self, chrom):
    p = random.random ()
    if p < self.mutation_probability:
        t = random.randint (1, self.chromosome_size)
        mask1 = 1 << (t - 1)
        mask2 = chrom & mask1
        if mask2 > 0:
            chrom = chrom & (~mask2)
        else:
            chrom = chrom ^ mask1
    return chrom
@

mutate 函式可以將一個染色體按照變異概率進行單點變異。例如:

\begin{verbatim}
1000 1101 1100 0010 0001 0110 1
\end{verbatim}

\noindent 在第 13 位發生變異,結果為:

\begin{verbatim}
1000 1101 1100 1010 0001 0110 1
\end{verbatim}

同交叉類似,也可以進行{\bf 多點變異}。

\section{進化}

<<進化機制>>=
def evolve (self):
    indvs = self.individuals
    new_indvs = self.new_individuals
    
    # 計算適應度及選擇概率
    self.evaluate ()
    
    # 進化操作
    i = 0
    while True:
        # 選擇兩個個體,進行交叉與變異,產生新的種群
        idv1 = self.select ()
        idv2 = self.select ()
        
        # 交叉
        (idv1_x, idv1_y) = (indvs[idv1][0], indvs[idv1][1])
        (idv2_x, idv2_y) = (indvs[idv2][0], indvs[idv2][1])
        (idv1_x, idv2_x) = self.cross (idv1_x, idv2_x)
        (idv1_y, idv2_y) = self.cross (idv1_y, idv2_y)
        
        # 變異
        (idv1_x, idv1_y) = (self.mutate (idv1_x), self.mutate (idv1_y))
        (idv2_x, idv2_y) = (self.mutate (idv2_x), self.mutate (idv2_y))
        
        (new_indvs[i][0], new_indvs[i][1])     = (idv1_x, idv1_y)
        (new_indvs[i+1][0], new_indvs[i+1][1]) = (idv2_x, idv2_y)
        
        # 判斷進化過程是否結束
        i = i + 2
        if i >= self.size:
            break
    
    # 更新換代
    for i in range (self.size):
        self.individuals[i][0] = self.new_individuals[i][0]
        self.individuals[i][1] = self.new_individuals[i][1]
    
@

[[evolve]] 函式可以實現種群的一代進化計算,計算過程分為三個步驟:

\begin{itemize}
\item 使用 [[evaluate]] 函式評估當前種群的適應度,並計算各個體的選擇概率。
\item 對於數量為 [[self.size]] 的 [[self.individuals]] 集合,迴圈 $\text{self.size}/ 2$ 次,每次從 [[self.individuals]] 中選出 2 個個體,對其進行交叉和變異操作,並將計算結果儲存於新的個體集合 [[self.new_individuals]] 中。
\item 用種群進化生成的新個體集合 [[self.new_individuals]] 替換當前個體集合。
\end{itemize}

如果迴圈呼叫 [[evolve]] 函式,那麼便可以產生一個種群進化的過程,如下:

<<進化機制>>=
def run (self):
    for i in range (self.generation_max):
        self.evolve ()
        print (i, max (self.fitness), sum (self.fitness)/self.size, 
               min (self.fitness))
@ 

[[run]] 函式根據種群最大進化世代數設定了一個迴圈。在迴圈過程中,呼叫 [[evolve]] 函式進行種群進化計算,並輸出種群的每一代的個體適應度最大值、平均值和最小值。

\section{開啟上帝模式}

下面的程式碼可以啟動種群進化過程:

<<啟動一個種群的進化過程>>=
if __name__ == '__main__':
    # 種群的個體數量為 50,染色體長度為 25,交叉概率為 0.8,變異概率為 0.1,進化最大世代數為 150
    pop = Population (50, 24, 0.8, 0.1, 150)
    pop.run ()
@

注意,因為個體交叉的需求,種群所包含的個體數量一般設為偶數。這個程式沒考慮個體數量為奇數的情況。

如果將以上所有出現的 python 程式碼依序組裝在一起,假設存為 test.py 檔案:

<<hello-ga.py>>=
import math, random

<<種群>>:
    <<種群初始化>>
    <<物競天擇機制>>
    <<染色體交叉機制>>
    <<染色體變異機制>>
    <<進化機制>>
    
<<啟動一個種群的進化過程>>
@ 

執行以下命令便可執行這個程式:

\begin{verbatim}
$ python3 test.py
\end{verbatim}

注意,這裡我們使用的是 python 3。如果你的系統中只安裝了 python 2,要讓程式能夠執行,需要在 test.py 的首行新增:

\begin{verbatim}
# -*- coding: utf-8 -*-
\end{verbatim}

然後將 [[evolve]] 函式中的 [[print]] 語句修改為:

\begin{verbatim}
print i, max (self.fitness), sum (self.fitness)/self.size, min (self.fitness)
\end{verbatim}

\section{結果}

如果使用命令:

\begin{verbatim}
$ python3 hello-ga.py > test.log
\end{verbatim}

那麼使用下面的 gnuplot 指令碼 test.gnu 可以繪製出種群的每一代最大適應度、平均適應度和最小適應度的變化情況。

\begin{verbatim}
#!/usr/bin/gnuplot
set term pngcairo
set size ratio 0.75
set output 'test.png'
plot "test.log" using 1:2 title "max" with lines, \
     "test.log" using 1:3 title "ave" with lines, \
     "test.log" using 1:4 title "min" with lines
\end{verbatim}

執行這個 gnuplot 指令碼,可以生成圖片檔案。

\begin{verbatim}
chmod +x ./test.gnu
./test.gnu
\end{verbatim}

圖 \ref{fig:第一次測試的結果} 中紅色的折線表示種群每一代個體中適應度最大值的變化情況,顯然,我們所得結果是比較接近 $f(x,\,y)$ 理論上的最大值 1.0。藍色折線反映了種群每一代最差個體適應度的變動情況,它的波動幅度看上去比較劇烈。如果將變異概率設為 0.4,那麼它看起來就會比較溫順一些,如圖 \ref{fig:第二次測試的結果} 所示。變異概率如果設定的越大,那麼藍色折線的波動幅度便會越小。圖 \ref{fig:第三次測試的結果} 顯示了比較極端的情況,此時變異概率設為 1.0。

在固定變異概率的條件下,可以用類似的方法觀察一下交叉概率對計算結果的影響。

\begin{figure}[H]
\centering
\includegraphics[width=8cm]{ga-test.png}
\caption{交叉概率為 0.8,變異概率為 0.1,種群的進化過程}
\label{fig:第一次測試的結果}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=8cm]{ga-test-1.png}
\caption{交叉概率為 0.8,變異概率為 0.4,種群的進化過程}
\label{fig:第二次測試的結果}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=8cm]{ga-test-2.png}
\caption{交叉概率為 0.8,變異概率為 1.0,種群的進化過程}
\label{fig:第三次測試的結果}
\end{figure}

\section{討論}

研究遺傳演算法的人證明了幾個定理。

\begin{theorem}
標準遺傳演算法不能收斂至全域性最優解。
\end{theorem}

本程式按照標準遺傳演算法實現的,從上面的幾幅圖也可以看出來,受交叉與變異的影響,種群的每一代個體的最大適應度都有可能在不斷變化。

\begin{theorem}
標準遺傳演算法,如果在選擇之前保留當前最佳個體,最終能收斂到全域性最優解。
\end{theorem}

對於本文所實現的遺傳演算法,只需要新增一個可以複製當前最佳個體資訊的函式,即可保證全域性最優解的收斂性,如下:

\begin{verbatim}
# 將該函式插入 Population 類中
    def reproduct_elitist (self):
        # 與當前種群進行適應度比較,更新最佳個體
        j = 0
        for i in range (self.size):
            if self.elitist['fitness'] < self.fitness[i]:
                j = i
                self.elitist['fitness'] = self.fitness[i]
        if (j > 0):
            self.elitist['chromosome'][0] = self.individuals[j][0]
            self.elitist['chromosome'][1] = self.individuals[j][1]
            self.elitist['age'] = self.age
\end{verbatim}

然後在 [[evlove]] 函式中呼叫 [[reporduct_elitist]] 函式:

\begin{verbatim}
# 修改後的 evolve 函式
    def evolve (self):
        indvs = self.individuals
        new_indvs = self.new_individuals
        
        # 計算適應度及選擇概率
        self.evaluate ()
        
        # 進化操作
        i = 0
        while True:
            # 選擇兩個個體,進行交叉與變異,產生新的種群
            idv1 = self.select ()
            idv2 = self.select ()
            
            # 交叉
            (idv1_x, idv1_y) = (indvs[idv1][0], indvs[idv1][1])
            (idv2_x, idv2_y) = (indvs[idv2][0], indvs[idv2][1])
            (idv1_x, idv2_x) = self.cross (idv1_x, idv2_x)
            (idv1_y, idv2_y) = self.cross (idv1_y, idv2_y)
            
            # 變異
            (idv1_x, idv1_y) = (self.mutate (idv1_x), self.mutate (idv1_y))
            (idv2_x, idv2_y) = (self.mutate (idv2_x), self.mutate (idv2_y))
            
            (new_indvs[i][0], new_indvs[i][1])     = (idv1_x, idv1_y)
            (new_indvs[i+1][0], new_indvs[i+1][1]) = (idv2_x, idv2_y)
            
            # 判斷進化過程是否結束
            i = i + 2
            if i >= self.size:
                break
        
        # 最佳個體保留
        self.reproduct_elitist ()
        
        # 更新換代
        for i in range (self.size):
            self.individuals[i][0] = self.new_individuals[i][0]
            self.individuals[i][1] = self.new_individuals[i][1]
\end{verbatim}

注意,我沒有將最佳個體儲存在種群的個體集合中,因為我覺得一個既不參與交叉也不參與變異的個體,是不能放在種群中的,它應當存放在歷史課本里。所以,我在 [[Population]] 類中設定了一個 [[elitist]] 的成員,用以記錄最佳個體對應的染色體、適應度及其出現的年代。當遺傳演算法結束後,這個最佳個體可作為目標函式的解。

\begin{theorem}
遺傳演算法所接受的引數有種群規模、適應度函式、染色體的表示、交叉概率、變異概率等,對於這些引數而言,不存在一個最佳組合,使它對於任何問題都能達到最優效能。
\end{theorem}

也就是說,成功的設計一個遺傳演算法的關鍵在於針對具體問題去選擇恰當的引數。如果不利用所求解一些問題的特定知識,那麼演算法的效能於我們採用何種引數沒有多大關係,情況可能會更糟糕。還要記住的是,對於單個問題,不存在最好的搜尋演算法。

 

新建一份文字檔案 hello-ga.nw,將上述程式碼複製到這份檔案中,然後執行以下命令

$ notangle -Rhello-ga.py hello-ga.nw > hello-ga.py

即可在當前目錄中生成 hello-ga.py 指令碼。也就是說,notangle 的功能是從文學程式設計的原始檔 hello-ga.nw 中提取 python 程式碼。

上述命令中的 -R 選項,是告訴 notangle 應該從一個叫做 hello-ga.py 的程式碼塊開始,而這個程式碼塊就是:

<<hello-ga.py>>= import math, random <<種群>>: <<種群初始化>> <<物競天擇機制>> <<染色體交叉機制>> <<染色體變異機制>> <<進化機制>> <<啟動一個種群的進化過程>>

要將 hello-ga.nw 檔案轉換為 LaTeX 文件,然後藉助 xelatex 生成 PDF 文件,可使用以下命令:

$ noweave -x hello-ga.nw | zhtex | > hello-ga.tex $ xelatex hello-ga

zhtex 是我自己寫的一份 Shell 指令碼,它藉助 sed 對 noweave 生成的 TeX 文件進行一些修改。這份指令碼的內容如下:

#!/bin/bash sed -e 's/\\documentclass{article}/\ \\documentclass[adobefonts]{ctexart}\ \\usepackage{amsmath}\ \\usepackage{graphicx, float}\ \\usepackage{cprotect}\ \\usepackage{booktabs}\ \\newtheorem{theorem}{定理}\ \\usepackage[top=2cm,bottom=2cm,left=2cm,right=2cm]{geometry}\n/g ; s/\\nwbegindocs{0}//g ; s/\\maketitle/\\maketitle\n\\nwbegindocs{0}/g'

生成的 hello-ga.pdf 文件,下載地址在:http://pan.baidu.com/s/1hq4QiGW

如果要用 noweb 來寫論文式程式,並希望能能生成含有數學公式、圖片、表格的 PDF 文件,你至少需要安裝 noweb 與 TeXLive。在 Gentoo 中,由於 TeXLive 整個包被拆分了,要讓 TeXLive 支援中文,需要安裝 texlive-langchinese。

noweb 中的 no,應該是作者的名字 Norman 的前兩個字母。

我為 noweb 的使用寫了一份簡單的指南,見『noweb 的用法』。

為什麼慢

如果你看了上文中的論文式程式設計程式碼以及所生成的 hello-ga.pdf 文件,可能會受到一些啟發,甚至在業餘時間裡也嘗試使用 noweb 來寫一些論文式程式。如果你真的這麼做了,很快就會發現,事情並沒有那麼美好。

因為在寫『論文』的過程中,你無法及時的驗證論文中的程式程式碼是否正確,只能等到論文式原始檔中包含了程式的一個完整的子集(即提取出的程式碼可被編譯或解釋執行),然後方有機會去測試程式程式碼的正確性。如果你堅持一次性的將論文寫完,然後再去驗證其中的程式程式碼是否正確,那時可能已經積攢了一大堆錯誤了。我們不是機器,我們大腦的容錯能力非常強,所以很多程式碼在我們看來是『正確』的,而編譯器或直譯器會很生氣的說 no!

我那個 hello-ga.nw 程式,其實是從幾年前我寫的一篇博文『移植』而成的,而那篇博文事實上是在我已經寫出來 hello-ga.py 指令碼之後才寫的。

我不認為真的有人能在論文裡將程式寫正確。實踐論文式程式設計,最可行的辦法是先寫引言部分,然後全力以赴的去寫程式、驗證程式的正確性,最後再將自己所寫的程式碼論文化。所以,這個過程總是要比非文學程式設計方式多耗費 1~3 倍的時間。其實,這個過程與嚴肅的程式設計過程並沒有什麼本質區別,我們先思考,然後寫程式碼,最後再寫程式文件。文學程式設計的真正意義在於,它強調了文件的重要性——文件一直是程式猿最不想寫的東西,並統一了文件與程式碼——至少在形式上是這樣。

相關文章