「集體智慧編程」之第五章：「求最優解」的算法

優化

顯然剛開始我並不太明白這一章要講什麼，因爲根據「優化」這個詞，我還以爲是對函數進行優化之類的。後來，我才明白，這一章在要講求最最優解的算法。由於我曾在老師的算法課上講過遺傳算法，遺傳算法就是用來求最優解的算法，所以我忽然明白了這個優化具體是指什麼。

制定旅行計劃的例子

情況描述

顯然上述內容都是廢話，對於這一章，關鍵還是要怎麼學好那些優化算法。我自己是通過舉例子才明白了本章的意圖，所以我突然覺得舉例子是一個很好的方式。所以我們來引出本章的第一個舉例：爲一次家族旅遊制定計劃。

當然本次旅遊也不是說你想怎麼樣就怎麼樣的，有很多限制條件，在這些限制條件之下，怎麼安排這次計劃，具體來說，就每個人出行的時間、乘坐航班的時間、是否轉機、租用車輛的時間。綜合來說，我們並不是要某一個人感到他自己方便了，而是要所有的人在相互磨合之下，所有人都方便了，並且要保證成本消耗最低，比如機票的價格不同吧？比如租車時間的長短，這是因爲租車要花錢吧。顯然這一段是一個抽象的描述，下面，我們來看看如何具體到底有什麼要求？注意，會非常具體，具體到每一個航班的始發時間和價格。

準備工作

首先，要明白是六個人（親戚吧），他們本來分散在美國的全國各地，然後呢？他們約好了某一天，一起到紐約旅遊，本地約定一天往返，他們就旅遊這麼一天。所以，我們可以寫出如下
代碼：

 #將要去旅行的人，第一個是名字，第二個是目前所在地
 people = [('Seymour','BOS'),
           ('Franny','DAL'),
           ('Zooey','CAK'),
           ('Walt','MIA'),
           ('Buddy','ORD'),
           ('Les','OMA')]
 #他們都要到美國的紐約集合，這是旅行的目的地
 #LGA是紐約的機場
 destination = 'LGA'

從他們所在地到紐約是乘坐飛機，但是飛機有不同的航班，航班的價格也不一樣。書中爲我們準備了一份航班列表：schedule.txt，供我們使用，我節選其中一部分進行講解,如下所示

LGA,OMA,18:25,20:34,205
OMA,LGA,18:12,20:17,242

他們用逗號分隔，依次是始發地、目的地、飛行時間、到達時間和價格。由於之前我們已經說了，這個旅遊就是一天的時間的，一天之內往返，所以航班信息只有時間，沒有日期。

我們必然會在代碼中會使用這些航班信息，因爲我們要針對每一個人，選出一對航班，讓他從當前位置飛到紐約，再從紐約飛回他的家。所以，先讓我們把這份文件讀進代碼，我們決定放在一個字典內，再次強調，字典就是鍵值對。
其中setdefault的函數，比較難理解，可以 參考這裏
代碼如下：

 #我們將所有的航班信息讀到一個字典內，以起止點爲鍵，其他的爲值
 flights={}


 for line in file('schedule.txt'):
     origin,dest,depart,arrive,price = line.strip().split(',')


     #其中setdefault是作爲字典類的一個方法，主要的作用是應對一種情況：當一個鍵對應多個值,每一個值是一個元組，多個元祖組成一個列表，
     #說白了，一個鍵對應一個列表，列表內有許多的元祖，一個元祖代表一個航班信息
     #setdefault的含義是："如果沒有，就設置"，如果有，就添加。
     flights.setdefault((origin,dest),[])
     #現在我們一行一行的把航班信息加入進去
     #我認爲，我們是對一個列表作爲值，所以最外面是一個[]
     #裏面跟了一個元祖，這個元祖就是鍵
     flights[(origin,dest)].append((depart,arrive,int(price)))
     #通過print flight，我們可以看出在該字典中具體的存儲情況。一個起止點作爲一個值，對應了一個列表，而列表中有許多元祖，每一個元祖都是一個航班信息的起飛時間、降落時間和價格

此時，我們需要加入一個函數，函數是將幾點幾分，轉換爲從凌晨0點開始經過了多少分鐘，需要這個函數的原因：幾點幾分的比較非常不便，而單純的分鐘數，使得飛機的飛行時間、候機時間變得容易計算。

代碼：

 def getminutes(t):
     #strptime是將特定的時間格式轉化爲元組。
     x=time.strptime(t,'%H:%M')
     return x[3]*60+x[4]
 #我覺得x[0]x[1]x[2]是年月日

下面我們要規定一種方式，用這種方式用來簡單的表達某一個人選擇了哪兩個航班。實際上這就是我們的解，我們求的最優解，就是用這種方式來表達。從書中讀到，我們要注意對類似的問題、類似的解做到通用，顯然這個思維非常重要，但是目前對我來說，缺少太多實踐經驗，所以不容易。在題本中，我們用一組數組數列來表示，數組在Python裏面叫做列表。如下所示:

[1,4,3,2,7,3,6,3,2,4,5,3]

我們有6個人，每人兩個航班，所有上面一個有12個數字，依次的第一個人的起飛航班和返回航班，而1代表，這個人將會乘坐這天中從某地到紐約的第二趟航班，因爲0代表了第一趟。爲什麼這種方式會便捷呢？因爲在flights字典中，我們已經把起止點相同的航班合併爲了一個列表，而列表當然可以用其下標0、1下表示。


顯然，如果我們每次都這麼讀這一串數字是非常辛苦的，所以我們決定寫一個函數，這個函數接受這一串數字列表，然後將人名、和起飛時間，起飛地點等信息都完整的打印出來。
代碼：

 def printschedule(r):
     for d in range(len(r)/2):#針對每一個人打印兩個航班信息，我們只用重複r的一半次
         name=people[d][0]
         origin= people[d][1]
         out=flights[(origin,destination)][r[2*d]]
         ret=flights[(destination,origin)][r[2*d+1]]
         print '%10s%10s %5s-%5s $%3s %5s-%5s $%3s' % (name,origin,out[0],out[1],out[2],ret[0],ret[1],ret[2])

我上面這段代碼改錯改了很久，主要就是少了符號"："還有就是少了'’的s。

成本函數

下面我們要做一個函數，叫做成本函數。顧名思義，就是我們做一件事的成本和代價。該函數有一個輸入：就是航班信息，有一個輸出，就是一個數字，數字越大代表代價越高。
書中提到：成本函數是用最優化算法的關鍵，這是因爲它往往難以確定，到底哪些是成本，成本的大小有多少？通過這個例子我相信大家有着更深刻的體會。下面請一個一個看，有哪些成本和這些成本如何用金錢來度量。

• 飛機的價格：6個人乘坐12次飛機的機票費用，顯然我們希望價格越低越好。價格的金錢度量很直接。
• 飛行時間：每個人在乘坐飛機時所花費的時間，我們當然希望乘坐飛機的時間越短越好。每分鐘1美元。（代碼中忽略了）
• 等待時間：在機場等待其他成員到達的時間，當然我們希望這個時間越短越好。每分鐘1美元。
• 出發時間：早晨太早的飛機將會產生額外的成本，該成本就是要求旅遊的人減少睡眠時間。我們希望儘可能的不要早，要合適。所以這個數字應該是越晚越好。（代碼中忽略了）
• 汽車租用時間：集體租用一輛車，那麼因爲只旅遊一天而已，所以租車的時間最好要控制在二十四小時之內。超過之後會產生額外的成本。超過二十四小時則產生50美元的罰款。

上述這是書中作者找到的關鍵的成本因素和金額度量。我們按照自己覺得合理的方式隨意改動。

代碼如下：

 def schedulecost(sol):
     totalprice=0
     latestarrival=0#最晚到底時間
     earliestdep=24*60#最早離開時間，現在24*60是最好的情況，等一下會根據實際飛機情況發生改變


     for d in range(len(sol)/2):
         #得到每一個人的兩次航班的價格並且加入到總價格中
         origin = people[d][1]
         outbound = flights[(origin,destination)][int(sol[2*d])]
         returnf = flights[(destination,origin)][int(sol[2*d+1])]


         #把錢加入到總價格裏面去
         totalprice+=outbound[2]
         totalprice+=returnf[2]


         #根據實際情況改變最晚到達時間和最早離開時間
         if latestarrival<getminutes(outbound[1]):latestarrival=getminutes(outbound[1])
         if earliestdep>getminutes(returnf[0]):earliestdep=getminutes(returnf[0])


     #每一個人必須在機場等待直到最後一個來了才能出發
     #每一個人必須在旅遊結束時，爲了最早離開的人能夠趕上飛機，而來到機場等候
     totalwait=0
     for d in range(len(sol)/2):
         origin = people[d][1]
         outbound = flights[(origin,destination)][int(sol[2*d])]
         returnf = flights[(destination,origin)][int(sol[2*d+1])]


         totalwait+=latestarrival-getminutes(outbound[1])
         totalwait+=earliestdep-getminutes(returnf[0])


     if latestarrival<earliestdep: totalprice+=50
     #如果這樣看，在機場每多等一分鐘就是1美元，而租車如果超過24個小時就罰款50美元
     return totalprice+totalwait

現在我成功打印出成本了，現在目標就是找出正確的數字序列，使成本最低。根據書中解釋，一共是12次航班，每次有10種不同的航班。那麼可以得到10的12次方個組合，就大約1000億。如果我們將這1000億個逐一比較，肯定能找到最佳答案。但是在計算機上耗費的時間很長。

隨機搜索

所以我們需要求最優解。在正在開始學習最優化算法之前，讓我們先來學一個叫隨機搜索的東西。
隨機搜索不是最優化算法，但是我們用它來評估其他算法的優劣。實際上，他就是隨機產生結果。但是我們可以重複隨機產生一堆結果...用來幹什麼我就不知道了。下面我們直接看代碼吧，關鍵的部分我會在相應出解釋。

代碼：

 #domain代表隨機產生的數字的個數和每一個數字的範圍，是一個列表，列表裏每個元素裏面是一個元組，每個元組有2個元素，一個是上限，一個是下限
 #costf就是成本函數，
 #每次隨機產生一組結果的時候，我們將會使用costf進行一下測試，看看效果如何
 def randomptimize(domain,costf):
     best=999999999
     bestr=None
     for i in range(10000):#我們打算隨機產生1000次結果，從這1000次結果中選擇一個最好的
         #很顯然randint是產生在一定範圍內的隨機數,顯然由於下一句右邊等號裏的for，將會產生一個循環
         r=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]
         cost=costf(r)


         #每次得到成本我們都判斷一次，如果更低，我們就置換
         if cost<best:
             best=cost
             bestr=r
     return bestr

        
其中對domain這個理解需要比較到位，如下產生，打印出來看看就明白是什麼了

 >>> domain=[(0,9)]*(len(people)*2)
 >>> print domain
 [(0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9),
  (0, 9), (0, 9)]

ok，在做了這麼多準備工作之後，我們終於可以開始學習真在的求最優解的算法了。

爬山法

核心思維

通過隨機搜索我們可以發現，每次尋找的解都是跳躍性的，前一個尋找的解和後一個尋找的解沒有任何聯繫。所以在尋找最優解的過程中，效率不過，我們現在的想法是：讓每一個解都對尋找最優解提供一點點幫助。爬山法的核心思維是：從一個隨機解開始，然後在其臨近的解集中尋找更好地題解，在本題中，就是更低的成本。核心思維有了，但是如何拿到臨近解確實一個問題，或者說，如何設置臨近解的產生，更加確切，在本題中，我們針對一個解，其中的每一個數字都加一或者減一，這樣我們就能使得某一個人坐的飛機能夠稍早一些或者稍晚一些。之所以加一或減一能夠決定一個人稍稍早一些或者晚一些，這是因爲：

經過查看flights{}字典發現，原來每一個鍵都是按照時間順序排列的

如下所示，我記得我的代碼並沒有刻意爲之排序，看樣子是默認排序的。

 >>>
 {('LGA', 'CAK'): [('6:58', '9:01', 238), ('8:19', '11:16', 122), ('9:58', '12:56', 249), ('10:32', '13:16', 139), ('12:01', '13:41', 267), ('13:37', '15:33', 142), ('15:50', '18:45', 243), ('16:33', '18:15', 253), ('18:17', '21:04', 259), ('19:46', '21:45', 214)], ('DAL', 'LGA'): [('6:12', '10:22', 230), ('7:53', '11:37', 433), ('9:08', '12:12', 364), ('10:30', '14:57', 290), ('12:19', '15:25', 342), ('13:54', '18:02', 294), ('15:44', '18:55', 382), ('16:52', '20:48', 448), ('18:26', '21:29', 464), ('20:07', '23:27', 473)], ('LGA', 'BOS'): [('6:39', '8:09', 86), ('8:23', '10:28', 149), ('9:58', '11:18', 130), ('10:33', '12:03', 74), ('12:08', '14:05', 142), ('13:39', '15:30', 74), ('15:25', '16:58', 62), ('17:03', '18:03', 103), ('18:24', '20:49', 124), ('19:58', '21:23', 142)], ('LGA', 'MIA'): [('6:33', '9:14', 172), ('8:23', '11:07', 143), ('9:25', '12:46', 295), ('11:08', '14:38', 262), ('12:37', '15:05', 170), ('14:08', '16:09', 232), ('15:23', '18:49', 150), ('16:50', '19:26', 304), ('18:07', '21:30', 355), ('20:27', '23:42', 169)], ('LGA', 'OMA'): [('6:19', '8:13', 239), ('8:04', '10:59', 136), ('9:31', '11:43', 210), ('11:07', '13:24', 171), ('12:31', '14:02', 234), ('14:05', '15:47', 226), ('15:07', '17:21', 129), ('16:35', '18:56', 144), ('18:25', '20:34', 205), ('20:05', '21:44', 172)], ('OMA', 'LGA'): [('6:11', '8:31', 249), ('7:39', '10:24', 219), ('9:15', '12:03', 99), ('11:08', '13:07', 175), ('12:18', '14:56', 172), ('13:37', '15:08', 250), ('15:03', '16:42', 135), ('16:51', '19:09', 147), ('18:12', '20:17', 242), ('20:05', '22:06', 261)], ('CAK', 'LGA'): [('6:08', '8:06', 224), ('8:27', '10:45', 139), ('9:15', '12:14', 247), ('10:53', '13:36', 189), ('12:08', '14:59', 149), ('13:40', '15:38', 137), ('15:23', '17:25', 232), ('17:08', '19:08', 262), ('18:35', '20:28', 204), ('20:30', '23:11', 114)], ('LGA', 'DAL'): [('6:09', '9:49', 414), ('7:57', '11:15', 347), ('9:49', '13:51', 229), ('10:51', '14:16', 256), ('12:20', '16:34', 500), ('14:20', '17:32', 332), ('15:49', '20:10', 497), ('17:14', '20:59', 277), ('18:44', '22:42', 351), ('19:57', '23:15', 512)], ('LGA', 'ORD'): [('6:03', '8:43', 219), ('7:50', '10:08', 164), ('9:11', '10:42', 172), ('10:33', '13:11', 132), ('12:08', '14:47', 231), ('14:19', '17:09', 190), ('15:04', '17:23', 189), ('17:06', '20:00', 95), ('18:33', '20:22', 143), ('19:32', '21:25', 160)], ('ORD', 'LGA'): [('6:05', '8:32', 174), ('8:25', '10:34', 157), ('9:42', '11:32', 169), ('11:01', '12:39', 260), ('12:44', '14:17', 134), ('14:22', '16:32', 126), ('15:58', '18:40', 173), ('16:43', '19:00', 246), ('18:48', '21:45', 246), ('19:50', '22:24', 269)], ('MIA', 'LGA'): [('6:25', '9:30', 335), ('7:34', '9:40', 324), ('9:15', '12:29', 225), ('11:28', '14:40', 248), ('12:05', '15:30', 330), ('14:01', '17:24', 338), ('15:34', '18:11', 326), ('17:07', '20:04', 291), ('18:23', '21:35', 134), ('19:53', '22:21', 173)], ('BOS', 'LGA'): [('6:17', '8:26', 89), ('8:04', '10:11', 95), ('9:45', '11:50', 172), ('11:16', '13:29', 83), ('12:34', '15:02', 109), ('13:40', '15:37', 138), ('15:27', '17:18', 151), ('17:11', '18:30', 108), ('18:34', '19:36', 136), ('20:17', '22:22', 102)]}
 >>>

代碼

核心思維結束了，讓我們直接來看代碼吧

 def hillclimb(domain,costf):
     #先創建一個隨機的解
     sol=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]
     while 1:#持續一個循環，直到在一次對每一個解減一或者加一之後沒有任何改變時，就break
         neighbors=[]
         for j in range(len(domain)):
             #解中的每一個元素都都會加一或者減一，加一產生一個解集，減一產生一個解集
             if sol[j]>domain[j][0]:
                 #如果很熟悉
                 neighbors.append(sol[0:j]+[sol[j]-1]+sol[j+1:])
             if sol[j]<domain[j][1]:
                 neighbors.append(sol[0:j]+[sol[j]+1]+sol[j+1:])
         #neighbors已經組裝好了，現在在裏面找最優解
         current=costf(sol)
         best=current
         for j in range(len(neighbors)):
             cost=costf(neighbors[j])
             if cost<best:
                 best=cost
                 sol=neighbors[j]
         if best == current:
             break
     return sol

可以調用如下語句直接產生一個由爬山法制作出來的結果，當然請保證其他函數的存在。

 domain=[(0,9)]*(len(people)*2)
 sol=hillclimb(domain,schedulecost)
 printschedule(sol)
 print schedulecost(sol)

難點

其中有一個難點：sol[0:j]+[sol[j]-1]+sol[j+1:]

我做了個測試，看測試的還很清楚就能明白這個什麼東西，加三個列表加起來，我覺得這個方法有點叼哇。最後一個錯誤請記住，要將三個列表相加的話，不要企圖用數字。

 >>> sol=[1,4,3,2,7,3,6,3,2,4,5,3]
 >>> a=[]
 >>> a.append(sol[0:0]+[sol[0]-1]+sol[1:])
 >>> print a
 [[0, 4, 3, 2, 7, 3, 6, 3, 2, 4, 5, 3]]
 >>> a.append(sol[0:0]+sol[0]-1+sol[1:])
 Traceback (most recent call last):
   File "<stdin>", line 1, in <module>
 TypeError: can only concatenate list (not "int") to list
 >>>

缺點

下面談一下爬山法的缺點。我認爲書本上的一幅圖，非常好體現了爬山法缺點。簡單來講，爬山法只能找到局部最優解，不能找到全局最優解，而我們要找的當然就是全局最優解。當然我們可以重複使用爬山法，都使用不同的初始狀態。

模擬退火算法

原理

受物理領域的一個啓發。說白了很簡單，根據剛剛所說的爬山法，我們很容易陷入局部最優解，這是因爲每次我們都只允許成本小的存在（成爲當前最優解，接着進行下一個循環）。我們沒有給算法一個跳出這個局部的機會，所以，我們這個模擬退火算法的關鍵就是：如果成更高，我們讓其有可能成爲當前的最優解。關鍵的關鍵就是在可能的程度上，我們要求：剛開始的時候可能大，隨着時間的推移（在模擬退火算法中是溫度的降低，當溫度降低到某個時刻），可能性越來越小，也就是隻接受成本更低的數列。這就有了跳出局部最優解的機會。那麼可能的程度，也就接受成本更大的概率，就使用一個公式來判斷：

爲什麼剛開始的概率大呢？因爲溫度高，所以指數接近了0，所以概率幾乎爲1,隨着溫度的減少，高成本與低成本之間的差值將會越來越重要，差異越大，概率越低。這注定了：該算法只會接受稍稍高一點成本的解，而不會接受成本高出許許多多的解。

代碼

具體執行代碼如下：

 def annealingoptimize(domain,costf,T=10000.0,cool=0.98,step=1):
     #和爬山法一樣，先產生一個隨機解，然後一切的改變都從這個隨機解開始
     vec=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]

     while T>0.5:
         #產生一個隨機數，決定這次改變是改變數列中的哪一個隨機數
         i=random.randint(0,len(domain)-1)


         #選擇一個改變的方向，也就是說是增加還是減少
         dir=random.randint(-step,step)


         #複製隨機解，然後對隨機解進行改變，然後判斷到底新的解好，還是後來產生的解好
         vecb=vec[:]
         vecb[i]+=dir
         #這一段主要還是不讓它超不過了最大最小值的限制
         if vecb[i]<domain[i][0]:vecb[i]=domain[i][0]
         elif vecb[i]>domain[i][1]:vecb[i]=domain[i][1]


         #計算新產生的兩次解的成本，然後對成本進行比較
         ea=costf(vec)
         eb=costf(vecb)

         #or後面：表示接受更差的結果。仔細想想，原來概率的表示是如此完成的，注意前一個random()產生的數是在0到1之間。
         if(eb<ea or random.random()<pow(math.e,-(eb-ea)/T)):
             vec=vecb


         #沒經過一次循環，改變溫度，溫度一改變，就會改變循環的次數和接受更差解的概率
         #按一定比例降溫
         T=T*cool


     return vec

使用一下代碼可以執行一次：

 domain=[(0,9)]*(len(people)*2)
 sol=annealingoptimize(domain,schedulecost)
 printschedule(sol)
 print schedulecost(sol)

經過我的驗證，確實可以產生更不錯的低成本，另外，偶爾也會得到一個較高成本的結果，記住一定要使用不同的參數多試一試。初始溫度，冷卻率，step的值等等。

遺傳算法

原理

遺傳算法是受生物學啓發而產生了的。該類算法首先產生一個組解，稱之爲種羣，對應種羣裏面的每一個解，我們都會計算其成本，顯然一組解，應了一組成本。我們將成本排序。

如下圖

依據這個排了序的結果，我們將會產生下一個子代。爲了使下一個子代的成本更低，我們首先選出目前這一代的優良品種，由於排序，很容易選擇前幾名，我們可以也約定好，選擇前多少名。

選了前幾名之後，顯然這一代的種羣數量還不夠，我們採用變異或者配對的方式來產生更多的子代，湊滿種羣的數量。

變異的做法就是指對一個解裏面某一個小部分進行小的改變，如下圖所示。6改爲5,0改成1，當然這肯定是有實際含義的對吧。

配對又稱交叉，我們將兩個優良品種各取一部分，組成一個新的解。當然，實際上，要根據實際情況的不同而進行調整。如下圖所示：

代碼

然後我們湊夠了下一代種羣，之後又排序，重複上述過程，每一重複一次，我們爲迭代一次，迭代的次數可以自己選定。
具體代碼如下：

 #popsize：一個種羣的規模大小
 #mutprob：種羣中進行變異，而不進行配對的概率。
 #elite：在當前種羣中被認爲優秀的子代的比例，優秀的子代可以直接傳入下一代
 #maxiter：迭代運行多少代


 def geneticoptimize(domain,costf,popsize=50,step=1,mutprob=0.2,elite=0.2,maxiter=100):
     #方法中還在定義方法
     #變異操作
     def mutate(vec):
         i=random.randint(0,len(domain)-1)
         #完成第增加或減少的概率各一半
         if random.random()<0.5 and vec[i]>domain[i][0]:
             return vec[0:i]+[vec[i]-step]+vec[i+1:]
         elif vec[i]<domain[i][1]:
             return vec[0:i]+[vec[i]+step]+vec[i+1:]
         else: return vec
     #交叉操作：貌似用python編程是好快的說，我感覺比較複雜的句子只要兩句麼，還是我c/c++沒學好
     def crossover(r1,r2):
         #爲什麼減2，其實想把這個一個數字列表劃分爲兩段，再各取一半
         i=random.randint(1,len(domain)-2)
         return r1[0:i]+r2[i:]

 #構造初始種羣
     pop=[]
     for i in range(popsize):
         vec = [random.randint(domain[i][0],domain[i][1])for i in range(len(domain))]
         pop.append(vec)

     #每一代有多少優勢物種，我們需要保留
     topelite=int(elite*popsize)
     #主循環
     for i in range(maxiter):
         #print pop #但是如果不加這句會使下一句出現一個bug，就是傳過去的v是None,但是我講pop全部打印出來的話，又沒有問題
         scores=[(costf(v),v) for v in pop]#列表裏面，每一個元素都是一個元組，每一個元組是由一個數字和一個列表構成
         scores.sort()
         ranked=[v for (s,v) in scores]



         #從中選擇我們覺得優勢的物種，然後保留
         pop=ranked[0:topelite]


         #如果種羣數量不夠，那麼我們使用變異或者配對，產生新的後代個體
         while len(pop)<popsize:
             #變異的概率，這是由我們設定的,雖然這裏是變異和配對只能選擇其一，但是我認爲是可以共同進行的
             if random.random()<mutprob:#如果這樣做，就是變異的少，交叉的多吧
                 #變異
                 c=random.randint(0,topelite)#注意是從優秀的子代中選出一個進行變異
                 pop.append(mutate(ranked[c]))
             else:
                 c1=random.randint(0,topelite)#從優秀的子代中選擇
                 c2=random.randint(0,topelite)#從優秀的子代中選擇
                 pop.append(crossover(ranked[c1],ranked[c2]))

         print scores[0][0]#注意打印的是成本

     return scores[0][1]#這裏返回的是航班序列

使用以下代碼可以打印個結果:

 domain=[(0,9)]*(len(people)*2)
 sol=geneticoptimize(domain,schedulecost)
 printschedule(sol)
 print schedulecost(sol)

總結

到了這裏算是一個階段性結束。我們主要對一個模型抽象出了數字列表，然後對 隨機搜索 、 模擬退火算法 、 遺傳算法 進行了初步的理解和實現。我們發現在這裏一系列問題中最困難的一步就是能否把問題潛在的解轉化爲數字列表。然後我個人覺得將整個問題轉換爲數字列表，或者能夠運用編程的方式來解決就是很困難的。


此外，上述的求最優解的方法能發揮的功效呢？和問題本身有着密切的聯繫。上述求最優解的方法都依賴於一個事實：最優解應該接近去其他優解。然而，顯示中可能存在這裏的問題，如下圖所示：


可以看出最優解的左右兩邊是否得陡峭，當產生了接近於最優解的解，我會認爲其不是優解而排除，最後我們只能陷入了局部最優解之中，上圖的左方。


如果放在航班例子中來看，就是我們如果從當天的第二次航班轉到第三次航班時，比將其轉到第八次航班更有可能降低成本。這是因爲航班有序排列，對每一個解集的一個航班加一或者減一效果肯定比減5、加6來的效果好。但是如果航班處於無序狀態，我們的求最優解的方式則不會有太大的作用，還不如隨機搜索。其中的航班處於無序狀態。什麼意思？起飛的航班必須有序呀。


書中下面一節將了通過真實的航班搜索來完成，這個單獨開一篇博客，因爲情況涉及調用外部api。會稍稍複雜一些。

有偏好條件的求最優解

書中叫做涉及偏好的優化。
抽象性的描述就是：如果將有限的資源分配給多個表達了偏好的人，並儘可能使他們滿意（或者根據他們的意願，儘可能的滿足他們的需要）。

這一次主要講另一個例子，這個例子是學生宿舍分配求最優解問題

但是該代碼可以很容易延伸到其他問題

• 在線紙牌遊戲中的玩家牌桌分配
• 大型編程項目中開發人員的bug分配
• 家庭成員的家務分配

核心思想 ：從個體中提取信息，將其組合產生最優解。

學生宿舍例子

具體描述

學生宿舍例子的具體描述是：一共有5間宿舍，每間宿舍有兩個隔間（顯然美國人的宿舍是一人一個隔間），由10名學生來競爭住所。每個人學生有一個第一志願和第二志願。我們先用代碼來描述一下基本情況，建一個py文件：dorm.py。代碼如下：

 import math
 import random
 import MyOptimization
 #一個代表宿舍集合的列表，注意每一個宿舍有兩個隔間
 dorms=['Zeus','Athena','Hercules','Bacchus','Pluto']

 #表述每個學生的第一志願以及第二志願
 prefs=[('Toby',('Bacchus','Hercules')),
         ('Steve',('Zeus','Pluto')),
         ('Andrea',('Athena','Zeus')),
         ('Sarah',('Zeus','Pluto')),
         ('Dave',('Athena','Bacchus')),
         ('Jeff',('Hercules','Pluto')),
         ('Fred',('Pluto','Athena')),
         ('Suzie',('Bacchus','Hercules')),
         ('Laura',('Bacchus','Hercules')),
         ('Neil',('Hercules','Athena'))]

我們不可能滿足所有人的需求。另外，該例子非常小巧，但是在實際生活中，也許人數非常多，而且每間宿舍的隔間可能有4個。

與航班問題相比，最大的困難在於題解了每間宿舍僅限兩個學生居住的約束條件。

解決思路

好吧，書上的內容我讀懂了，但是我現在未必有很深的領悟。對這個問題的解決方法就是：隨意產生一個解，但是這個解必定是有效的，滿足條件的。解雖然也是用一個數字列表來表示，但是我們卻不能直接看出某個學生選了什麼。但是我麼保證了這個數字列表在我們的規則之下確實能打印出獨一無二的結果。此外，我們每次用這個有效解計算成本，選擇成本最小的。核心思路如下：


將所有的宿舍用槽來表示，每個宿舍兩個槽。共有10個槽，相同宿舍的槽用同樣的數字表示，當某一個槽被選了之後，我們將會將其刪除出槽列表。如下所示：

[0,0,1,1,2,2,3,3,4,4]
如果第一個學生選了第一個，槽列表將剩餘：
[0,1,1,2,2,3,3,4,4]
依次，槽列表越來越少。

現在，我們需要一個domain，就是範圍，有這個範圍內產生的一系列數字，並不是代表學生選擇了哪一個宿舍，只是代表了學生選擇剩餘槽內的第幾號元素。比如：[0,0,0,0,0,0,0,0,0,0],它可以存在，爲什麼呢？因爲它只是代表了每一位學生都選擇了槽內的第一個元素。而是每次剩餘槽內的第一個元素並都代表第一個宿舍。
比如[1,1,1,1,1,1,1,1,1,1]確實無效的，爲什麼？因爲輪到最後一個選的時候，槽內只有一個了，他卻要選第二個，顯然已經數組越界了。


所以產生的這個domain的範圍確是：

 #[(0,9),(0,8),(0,7),(0,8).....(0,0)]
 domain=[(0,(len(dorms)*2)-i-1) for i in range(0,len(dorms)*2)]

代碼

下面就是打印出結果了，如果讀懂了上述的思維，代碼將會顯得很容易理解了：

 #最關鍵的就是傳入的vec這個序列的含義了
 #比如[0,0,0,0,0,0,0,0,0,0]
 #這並不代表某一個學生要選擇dorms裏面的第一個宿舍
 #它代表的含義是選擇了剩餘槽列表裏面的第一個。
 def printsolution(vec):
     slots=[]
     #每個宿舍兩個槽
     for i in range(len(dorms)):slots+=[i,i]
     print slots
     #遍歷每一個學生的安置情況
     for i in range(len(vec)):
         x = int(vec[i])
         print slots
         #從剩餘槽中選擇
         dorm=dorms[slots[x]]
         print prefs[i][0],dorm


         #刪除該槽
         del slots[x]

上述函數可以幫助我們打印出最終學生的選擇，下面一個成本函數。成本函數的參數和printsolution是一I樣的，我們針對每一個有效解（千萬注意：該有效解不能直接看出某個選擇選擇了什麼），計算了成本，如果成本低，我們就會保留。
成本函數：

 def dormcost(vec):
     cost=0
     #建立一個槽序列
     slots=[]
     #每個宿舍兩個槽
     for i in range(len(dorms)):slots+=[i,i]


     #遍歷每一個學生
     for i in range(len(vec)):
         x=int(vec[i])
         dorm=dorms[slot[x]]
         pref=prefs[i][1]
         #首選成本值爲0,次選成本值爲1
         if pref[0]==dorm:cost+=0
         elif pref[1]==dorm:cost+=1
         else cost+=3#不在選擇之列則成本爲3


         #刪除選中的槽
         del slots[x]


     return cost

有了上面兩個函數，我們只需要像之前解決航班問題一樣，將domain和成本函數傳入到求最優解的函數或者是隨機搜索就可以了。


經過我實際運行，確實產生了最優解，非常棒。
在文檔開頭加入：import MyOptimization


代碼：

 #[(0,9),(0,8),(0,7),(0,8).....(0,0)]
 domain=[(0,(len(dorms)*2)-i-1) for i in range(0,len(dorms)*2)]
 sol=MyOptimization.geneticoptimize(domain,dormcost)
 printsolution(sol)

知識遷移能力很重要，希望這道題的解答過程能夠深深地印在我的腦海了。

網絡可視化

含義

一切高深的問題用例子來講都可以講的很簡單，這是一個在社交網絡的運用，我們想搞清楚人們之間的聯繫，如果用圖可以直接看出來將是非常方便的。如下圖所示：

這是一個有點亂的網格，雖然能看出來誰是誰的朋友，但是顯然看不出一些關鍵人物，對比另外一幅圖。

ok，這一小節，主要就是將怎麼利用數據做出上一幅圖。當然其中會使用到求最優解的算法。爲什麼？因爲我們爲了看清楚一幅圖，我們總希望交叉線最少，我們總希望角度最好大一點...等等，所有在可以產生的那麼多圖中，我們想找到我們能看的最清楚的那一幅圖

準備工作

我們需要一點數據，來完成我們這次的學習，如下圖所示：
import math
people=['Charlie','Augustus','Veruca','Violet','Mike','Joe','Willy','Miranda']
links=[('Augustus','Willy')
       ('Mike','Joe')
       ('Miranda','Mike')
       ('Violet','Augustus')
       ('Miranda','Willy')
       ('Charlie','Mike')
       ('Veruca','Joe')
       ('Miranda','Augustus')
       ('Willy','Augustus')
       ('Joe','Charlie')
       ('Veruca','Augustus')
       ('Miranda','Joe')
    ]


我們要做的事就是：開發一個程序，能夠讀出上面的數據，主要是誰是誰的好朋友，然後畫出一個視覺上方便的圖。

核心思維

畫圖的核心思維或者算法就是 質點彈簧法（mass-and-spring algorithm）.算法受物理學啓發:各個節點相互之間都有推力，將彼此分離，但是如果彼此有聯繫，則拉近一點點。這樣的話，沒有關聯的點就會被推在外面，而彼此連接緊密的節點在拉在了一起。
然而，質點彈簧法無法避免交叉線，交叉線越多，觀察越困難。然而，這就是我們需要求最優解函數的地方。我們可以設置一個成本函數。最簡單的辦法：成本函數就是計算交叉線的個數。

如何計算交叉線的個數呢?有兩個關鍵關鍵1就是如何表達題解呢？我們用座標的方式的，每一個人都有一個x座標和一個y座標，所以把這些所有的座標一次放在一個數字列表中就可以，比如：
sol=[120,200,250,125.....
這表示第一個人charlie位置座標（120,200）的地方。
關鍵2就是計算交叉點的方式使用了一個公式：分數值。兩條線，不平行肯定交叉，但是我們題中出現的線段，所以我們主要考察線段是否交叉。對兩個線段分別求出分數值：如果分數值的範圍在0到1之間，那麼就交叉，不然就不交叉。
成本函數的接受一個題解，然後遍歷所有的一對線段，判斷是否交叉，每交叉一次，我們就交成本增加1.

請注意，在這裏我們使用求最優解算法的地方就在於求出題解（各個人的座標），而每一個題解都有成本，成本就是交叉數，我們要請出成本最低的題解，

代碼

接着來讓我們看一下求成本的代碼（也就是交叉線）：

 def crosscount(v):
     #將數字序列轉化爲一個person:(x,y)的字典

     loc=dict([(people[i],(v[i*2],v[i*2+1])) for i in range(0,len(people))])
     total=0



     #遍歷每一對連線
     for i in range(len(links)):
         for j in range(i+1,len(links)):


             #獲得座標位置
             (x1,y1),(x2,y2)=loc[links[i][0]],loc[links[i][1]]
             (x3,y3),(x4,y4)=loc[links[j][0]],loc[links[j][1]]

             #下列幾個算式涉及一個數學公式。不要太在意其複雜的表達方式，它只是遵循了數學公式
             den=(y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)


             #如果兩線平行，則den==0
             if den==0:continue


             #否則，ua與ub就是兩條交叉線的分數值
             #如果要顯示成小數點那樣的形式，必須要用float，這個很重要。這是與書上不同的地方
             ua=float(((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3)))/float(den)
             ub=float(((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3)))/float(den)

             #如果兩條線的分數值介於0和1之間，則兩線彼此交叉
             if ua>0 and ua<1 and ub>0 and ub<1:
                 total+=1
     return total

上述與書中不同的主要就是那個float，不加float的話，產生的數是四捨五入之後的數。

接着，我們可以用如下代碼產生一個優秀的題解：
注意引入MyOptimization.py

 #因爲我們決定將我們的座標系建立在一個400*400的像素圖中，下面寫10-370是爲了留點邊緣，
 #所以我們可以可以知道每一個座標的x或者y的範圍都是10到370，而每一個人有x和y座標，所以要乘以2
 domain=[(10,370)]*(len(people)*2)
 sol=MyOptimization.geneticoptimize(domain,crosscount)
 print sol

我們現在已經求出最優解了，我們還想把它畫出來。


畫出來必須使用一個庫文件，說白了就畫圖的唄。庫的名字叫做 Python Imaging Library 。安裝過程很簡單，直接去 下載 ，然後對應自己python版本下載一個，當然也要注意自己的平臺。
然後就是代碼，其實代碼部分非常簡單（不過讓自己寫還真的未必寫的出來）

 from PIL import Image,ImageDraw
 def drawnetwork(sol):
     #建立image對象
     img=Image.new('RGB',(400,400),(255,255,255))
     draw=ImageDraw.Draw(img)


     #建立一個字典，準備打印位置，和上面計算成本一個意思
     pos=dict([(people[i],(sol[i*2],sol[i*2+1])) for i in range(0,len(people))])


     #劃線，也就是說不會刻意去畫點，全是線。
     for (a,b) in links:
         draw.line((pos[a],pos[b]),fill=(255,0,0))


     #把人物名字繪製出來，看着人名就知道哪裏是個點了吧
     for n,p in pos.items():
         draw.text(p,n,(0,0,0))
     img.show()

由於每次產生的題解都是不同的，所以每次畫出來的圖都不一樣，但是現在可以看出來確實非常有趣了。下面貼一幅產生的圖


但是這個圖很醜（實際上很不錯，我產生了很多次才產生出這個圖的），關鍵是

• 兩線夾角非常小
• 兩點間的間距很近

爲了畫出好看的圖，我們想要杜絕上述兩種情況，以一個爲例，如果我們不想要兩個點靠的太近，那麼我們可以裏面的每對點進行遍歷，如果像素靠的太近了，我們就加一點成本。在本題中，如果兩點靠近的50個像素。就多加點成本，具體加了多少請看代碼。

代碼是在crosscount函數中先添加了一部分（在最後return之前）：

 def crosscount(v):
     #將數字序列轉化爲一個person:(x,y)的字典

     loc=dict([(people[i],(v[i*2],v[i*2+1])) for i in range(0,len(people))])
     total=0



     #遍歷每一對連線
     for i in range(len(links)):
         for j in range(i+1,len(links)):


             #獲得座標位置
             (x1,y1),(x2,y2)=loc[links[i][0]],loc[links[i][1]]
             (x3,y3),(x4,y4)=loc[links[j][0]],loc[links[j][1]]

             #下列幾個算式涉及一個數學公式。不要太在意其複雜的表達方式，它只是遵循了數學公式
             den=(y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)


             #如果兩線平行，則den==0
             if den==0:continue


             #否則，ua與ub就是兩條交叉線的分數值
             #如果要顯示成小數點那樣的形式，必須要用float，這個很重要。這是與書上不同的地方
             ua=float(((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3)))/float(den)
             ub=float(((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3)))/float(den)

             #如果兩條線的分數值介於0和1之間，則兩線彼此交叉
             if ua>0 and ua<1 and ub>0 and ub<1:
                 total+=2
     for i in range(len(people)):
         for j in range(i+1,len(people)):
             #獲得兩個結點的位置
             (x1,y1),(x2,y2)=loc[people[i]],loc[people[j]]


             #計算兩結點之間的距離
             dist=math.sqrt(math.pow(x1-x2,2)+math.pow(y1-y2,2))
             #對間距小於50個像素的結點進去判罰。
             if dist<150:total+=(1.0-(dist/150.0))
     return total

然後我也改了改其中的數據，產生了兩幅還不錯的圖。

總結

算法的應用場合也許有很多，但是看你是否去用心思考。

所有的難點在於：

• 確定題解的表示法
• 成本函數的確定

實際上，我們在處理三個類型的例子的時候，優化算法其實只寫了一個。但是針對不同的問題，費了很多周折才完成了成本函數和題解的確定。

具體來說，什麼樣的問題能夠使用求最優解的算法來解決呢？書中提出兩點

• 問題本身有一個定義好的成本函數
• 相似的解會產生相似的結果

雖然不是所有此類問題都可以用求最優解的算法來解決，但是我們都可以試着研究一下。

對項目的思考

我的項目，哪裏能用最優化算法？如果我能想到，我覺得，我很厲害。哈哈

不過暫時沒想的很完整。

會不會當我把幾個推薦引擎設計好之後，針對某一個用戶，求出不同推薦引擎所佔的比重，這個時候會不會用一下求最優解的算法呢？

那恐怕要先回答成本函數怎麼定？

這個感覺得依據大量用戶的操作數據吧。我的意思收集了大量一個用戶操作數據之後，就可以做。因爲比如引擎A產生的歌曲，用戶點了刪除，那麼說明了....什麼？

好像有點複雜。再想想吧。

全部源代碼

一個例子一個，共三份

旅遊例子

 # -*- coding: cp936 -*-
 import time
 import random
 import math
 #將要去旅行的人，第一個是名字，第二個是目前所在地
 people = [('Seymour','BOS'),
           ('Franny','DAL'),
           ('Zooey','CAK'),
           ('Walt','MIA'),
           ('Buddy','ORD'),
           ('Les','OMA')]
 #他們都要到美國的紐約集合，這是旅行的目的地
 #LGA是紐約的機場
 destination = 'LGA'

 #我們將所有的航班信息讀到一個字典內，以起止點爲鍵，其他的爲值
 flights={}

 for line in file('schedule.txt'):
     origin,dest,depart,arrive,price = line.strip().split(',')

     #其中setdefault是作爲字典類的一個方法，主要的作用是應對一種情況：當一個鍵對應多個值,每一個值是一個元組，多個元組組成一個列表，
     #說白了，一個鍵對應一個列表，列表內有許多的元組，一個元組代表一個航班信息
     #setdefault的含義是："如果沒有，就設置"，如果有，就添加。
     flights.setdefault((origin,dest),[])
     #現在我們一行一行的把航班信息加入進去
     #我認爲，我們是對一個列表作爲值，所以最外面是一個[]
     #裏面跟了一個元組，這個元組就是鍵
     flights[(origin,dest)].append((depart,arrive,int(price)))
     #通過print flight，我們可以看出在該字典中具體的存儲情況。一個起止點作爲一個值，對應了一個列表，而列表中有許多元組，每一個元組都是一個航班信息的起飛時間、降落時間和價格



 def getminutes(t):
     #strptime是將特定的時間格式轉化爲元組。
     x=time.strptime(t,'%H:%M')
     return x[3]*60+x[4]
 #我覺得x[0]x[1]x[2]是年月

 def printschedule(r):
     for d in range(len(r)/2):#針對每一個人打印兩個航班信息，我們只用重複r的一半次
         name=people[d][0]
         origin= people[d][1]
         out=flights[(origin,destination)][r[2*d]]
         ret=flights[(destination,origin)][r[2*d+1]]
         print '%10s%10s %5s-%5s $%3s %5s-%5s $%3s' % (name,origin,out[0],out[1],out[2],ret[0],ret[1],ret[2])

 def schedulecost(sol):
     totalprice=0
     latestarrival=0#最晚到底時間
     earliestdep=24*60#最早離開時間，現在24*60是最好的情況，等一下會根據實際飛機情況發生改變
     for d in range(len(sol)/2):
         #得到每一個人的兩次航班的價格並且加入到總價格中
         origin = people[d][1]
         outbound = flights[(origin,destination)][int(sol[2*d])]
         returnf = flights[(destination,origin)][int(sol[2*d+1])]

         #把錢加入到總價格裏面去
         totalprice+=outbound[2]
         totalprice+=returnf[2]

         #根據實際情況改變最晚到達時間和最早離開時間
         if latestarrival<getminutes(outbound[1]):latestarrival=getminutes(outbound[1])
         if earliestdep>getminutes(returnf[0]):earliestdep=getminutes(returnf[0])

     #每一個人必須在機場等待直到最後一個來了才能出發
     #每一個人必須在旅遊結束時，爲了最早離開的人能夠趕上飛機，而來到機場等候
     totalwait=0
     for d in range(len(sol)/2):
         origin = people[d][1]
         outbound = flights[(origin,destination)][int(sol[2*d])]
         returnf = flights[(destination,origin)][int(sol[2*d+1])]

         totalwait+=latestarrival-getminutes(outbound[1])
         totalwait+=getminutes(returnf[0])-earliestdep

     if latestarrival<earliestdep: totalprice+=50
     #如果這樣看，在機場每多等一分鐘就是1美元，而租車如果超過24個小時就罰款50美元
     return totalprice+totalwait

 #domain代表隨機產生的數字的個數和每一個數字的範圍，是一個列表，列表裏每個元素裏面是一個元組，每個元組有2個元素，一個是上限，一個是下限
 #costf就是成本函數，
 #每次隨機產生一組結果的時候，我們將會使用costf進行一下測試，看看效果如何
 def randomptimize(domain,costf):
     best=999999999
     bestr=None
     for i in range(10000):#我們打算隨機產生1000次結果，從這1000次結果中選擇一個最好的
         #很顯然randint是產生在一定範圍內的隨機數,顯然由於下一句右邊等號裏的for，將會產生一個循環
         r=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]
         cost=costf(r)

         #每次得到成本我們都判斷一次，如果更低，我們就置換
         if cost<best:
             best=cost
             bestr=r
     return bestr

 def hillclimb(domain,costf):
     #先創建一個隨機的解
     sol=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]
     while 1:#持續一個循環，直到在一次對每一個解減一或者加一之後沒有任何改變時，就break
         neighbors=[]
         for j in range(len(domain)):
             #解中的每一個元素都都會加一或者減一，加一產生一個解集，減一產生一個解集
             if sol[j]>domain[j][0]:
                 #如果很熟悉
                 neighbors.append(sol[0:j]+[sol[j]-1]+sol[j+1:])
             if sol[j]<domain[j][1]:
                 neighbors.append(sol[0:j]+[sol[j]+1]+sol[j+1:])
         #neighbors已經組裝好了，現在在裏面找最優解
         current=costf(sol)
         best=current
         for j in range(len(neighbors)):
             cost=costf(neighbors[j])
             if cost<best:
                 best=cost
                 sol=neighbors[j]
         if best == current:
             break
     return sol

 def annealingoptimize(domain,costf,T=10000.0,cool=0.98,step=1):
     #和爬山法一樣，先產生一個隨機解，然後一切的改變都從這個隨機解開始
     vec=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]

     while T>0.5:
         #產生一個隨機數，決定這次改變是改變數列中的哪一個隨機數
         i=random.randint(0,len(domain)-1)

         #選擇一個改變的方向，也就是說是增加還是減少
         dir=random.randint(-step,step)

         #複製隨機解，然後對隨機解進行改變，然後判斷到底新的解好，還是後來產生的解好
         vecb=vec[:]
         vecb[i]+=dir
         #這一段主要還是不讓它超不過了最大最小值的限制
         if vecb[i]<domain[i][0]:vecb[i]=domain[i][0]
         elif vecb[i]>domain[i][1]:vecb[i]=domain[i][1]

         #計算新產生的兩次解的成本，然後對成本進行比較
         ea=costf(vec)
         eb=costf(vecb)

         #or後面：表示接受更差的結果。仔細想想，原來概率的表示是如此完成的，注意前一個random()產生的數是在0到1之間。
         if(eb<ea or random.random()<pow(math.e,-(eb-ea)/T)):
             vec=vecb

         #沒經過一次循環，改變溫度，溫度一改變，就會改變循環的次數和接受更差解的概率
         #按一定比例降溫
         T=T*cool

     return vec

 #popsize：一個種羣的規模大小
 #mutprob：種羣中進行變異，而不進行配對的概率。
 #elite：在當前種羣中被認爲優秀的子代的比例，優秀的子代可以直接傳入下一代
 #maxiter：迭代運行多少代

 def geneticoptimize(domain,costf,popsize=50,step=1,mutprob=0.2,elite=0.2,maxiter=100):
     #方法中還在定義方法
     #變異操作
     def mutate(vec):
         i=random.randint(0,len(domain)-1)
         #完成第增加或減少的概率各一半
         if random.random()<0.5 and vec[i]>domain[i][0]:
             return vec[0:i]+[vec[i]-step]+vec[i+1:]
         elif vec[i]<domain[i][1]:
             return vec[0:i]+[vec[i]+step]+vec[i+1:]
         else: return vec
     #交叉操作：貌似用python編程是好快的說，我感覺比較複雜的句子只要兩句麼，還是我c/c++沒學好
     def crossover(r1,r2):
         #爲什麼減2，其實想把這個一個數字列表劃分爲兩段，再各取一半
         i=random.randint(1,len(domain)-2)
         return r1[0:i]+r2[i:]

 #構造初始種羣
     pop=[]
     for i in range(popsize):
         vec = [random.randint(domain[i][0],domain[i][1])for i in range(len(domain))]
         pop.append(vec)

     #每一代有多少優勢物種，我們需要保留
     topelite=int(elite*popsize)
     #主循環
     for i in range(maxiter):
         #print pop #但是如果不加這句會使下一句出現一個bug，就是傳過去的v是None,但是我講pop全部打印出來的話，又沒有問題
         scores=[(costf(v),v) for v in pop]#列表裏面，每一個元素都是一個元組，每一個元組是由一個數字和一個列表構成
         scores.sort()
         ranked=[v for (s,v) in scores]


         #從中選擇我們覺得優勢的物種，然後保留
         pop=ranked[0:topelite]

         #如果種羣數量不夠，那麼我們使用變異或者配對，產生新的後代個體
         while len(pop)<popsize:
             #變異的概率，這是由我們設定的,雖然這裏是變異和配對只能選擇其一，但是我認爲是可以共同進行的
             if random.random()<mutprob:#如果這樣做，就是變異的少，交叉的多吧
                 #變異
                 c=random.randint(0,topelite)#注意是從優秀的子代中選出一個進行變異
                 pop.append(mutate(ranked[c]))
             else:
                 c1=random.randint(0,topelite)#從優秀的子代中選擇
                 c2=random.randint(0,topelite)#從優秀的子代中選擇
                 pop.append(crossover(ranked[c1],ranked[c2]))

         print scores[0][0]#注意打印的是成本

     return scores[0][1]#這裏返回的是航班序列


 domain=[(0,9)]*(len(people)*2)
 #sol=geneticoptimize(domain,schedulecost)
 #printschedule(sol)
 #print schedulecost(sol)

宿舍分配例子

 # -*- coding: cp936 -*-
 import math
 import random
 import MyOptimization
 #一個代表宿舍集合的列表，注意每一個宿舍有兩個隔間
 dorms=['Zeus','Athena','Hercules','Bacchus','Pluto']

 #表述每個學生的第一志願以及第二志願
 prefs=[('Toby',('Bacchus','Hercules')),
         ('Steve',('Zeus','Pluto')),
         ('Andrea',('Athena','Zeus')),
         ('Sarah',('Zeus','Pluto')),
         ('Dave',('Athena','Bacchus')),
         ('Jeff',('Hercules','Pluto')),
         ('Fred',('Pluto','Athena')),
         ('Suzie',('Bacchus','Hercules')),
         ('Laura',('Bacchus','Hercules')),
         ('Neil',('Hercules','Athena'))]

 #最關鍵的就是傳入的vec這個序列的含義了
 #比如[0,0,0,0,0,0,0,0,0,0]
 #這並不代表某一個學生要選擇dorms裏面的第一個宿舍
 #它代表的含義是選擇了剩餘槽列表裏面的第一個。
 def printsolution(vec):
     slots=[]
     #每個宿舍兩個槽
     for i in range(len(dorms)):slots+=[i,i]
     #遍歷每一個學生的安置情況
     for i in range(len(vec)):
         x = int(vec[i])
         #從剩餘槽中選擇
         dorm=dorms[slots[x]]
         print prefs[i][0],dorm

         #刪除該槽
         del slots[x]



 def dormcost(vec):
     cost=0
     #建立一個槽序列
     slots=[]
     #每個宿舍兩個槽
     for i in range(len(dorms)):slots+=[i,i]

     #遍歷每一個學生
     for i in range(len(vec)):
         x=int(vec[i])
         dorm=dorms[slots[x]]
         pref=prefs[i][1]
         #首選成本值爲0,次選成本值爲1
         if pref[0]==dorm:cost+=0
         elif pref[1]==dorm:cost+=1
         else: cost+=3#不在選擇之列則成本爲3
         #刪除選中的槽
         del slots[x]
     print cost
     return cost

 #[(0,9),(0,8),(0,7),(0,8).....(0,0)]
 domain=[(0,(len(dorms)*2)-i-1) for i in range(0,len(dorms)*2)]
 sol=MyOptimization.geneticoptimize(domain,dormcost)
 printsolution(sol)

社交網絡繪圖

 # -*- coding: cp936 -*-
 import math
 import MyOptimization
 from PIL import Image,ImageDraw

 people=['Charlie','Augustus','Veruca','Violet','Mike','Joe','Willy','Miranda']
 links=[('Augustus','Willy'),
        ('Mike','Joe'),
        ('Miranda','Mike'),
        ('Violet','Augustus'),
        ('Miranda','Willy'),
        ('Charlie','Mike'),
        ('Veruca','Joe'),
        ('Miranda','Augustus'),
        ('Willy','Augustus'),
        ('Joe','Charlie'),
        ('Veruca','Augustus'),
        ('Miranda','Joe')]

 def crosscount(v):
     #將數字序列轉化爲一個person:(x,y)的字典

     loc=dict([(people[i],(v[i*2],v[i*2+1])) for i in range(0,len(people))])
     total=0


     #遍歷每一對連線
     for i in range(len(links)):
         for j in range(i+1,len(links)):

             #獲得座標位置
             (x1,y1),(x2,y2)=loc[links[i][0]],loc[links[i][1]]
             (x3,y3),(x4,y4)=loc[links[j][0]],loc[links[j][1]]

             #下列幾個算式涉及一個數學公式。不要太在意其複雜的表達方式，它只是遵循了數學公式
             den=(y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)



             #如果兩線平行，則den==0
             if den==0:continue

             #否則，ua與ub就是兩條交叉線的分數值
             #如果要顯示成小數點那樣的形式，必須要用float，這個很重要。這是與書上不同的地方
             ua=float(((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3)))/float(den)
             ub=float(((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3)))/float(den)

             #如果兩條線的分數值介於0和1之間，則兩線彼此交叉
             if ua>0 and ua<1 and ub>0 and ub<1:
                 total+=2
     for#否則，ua與ub就是兩條交叉線的分數值
             #如果要顯示成小數點那樣的形式，必須要用float，這個很重要。這是與書上不同的地方
             ua=float(((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3)))/float(den)
             ub=float(((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3)))/float(den)

             #如果兩條線的分數值介於0和1之間，則兩線彼此交叉
             if ua>0 and ua<1 and ub>0 and ub<1:
                 total+=2
     for i in range(len(people)):
         for j in range(i+1,len(people)):
             #獲得兩個結點的位置 i in range(len(people)):