“集体智慧编程”之第五章：“求最优解”的算法

优化

显然刚开始我并不太明白这一章要讲什么，因为根据“优化”这个词，我还以为是对函数进行优化之类的。后来，我才明白，这一章在要讲求最最优解的算法。由于我曾在老师的算法课上讲过遗传算法，遗传算法就是用来求最优解的算法，所以我忽然明白了这个优化具体是指什么。

制定旅行计划的例子

情况描述

显然上述内容都是废话，对于这一章，关键还是要怎么学好那些优化算法。我自己是通过举例子才明白了本章的意图，所以我突然觉得举例子是一个很好的方式。所以我们来引出本章的第一个举例：为一次家族旅游制定计划。

当然本次旅游也不是说你想怎么样就怎么样的，有很多限制条件，在这些限制条件之下，怎么安排这次计划，具体来说，就每个人出行的时间、乘坐航班的时间、是否转机、租用车辆的时间。综合来说，我们并不是要某一个人感到他自己方便了，而是要所有的人在相互磨合之下，所有人都方便了，并且要保证成本消耗最低，比如机票的价格不同吧？比如租车时间的长短，这是因为租车要花钱吧。显然这一段是一个抽象的描述，下面，我们来看看如何具体到底有什么要求？注意，会非常具体，具体到每一个航班的始发时间和价格。

准备工作

首先，要明白是六个人（亲戚吧），他们本来分散在美国的全国各地，然后呢？他们约好了某一天，一起到纽约旅游，本地约定一天往返，他们就旅游这么一天。所以，我们可以写出如下
代码：

[python] view plain copy

1. #将要去旅行的人，第一个是名字，第二个是目前所在地  
2. people = [('Seymour','BOS'),  
3.           ('Franny','DAL'),  
4.           ('Zooey','CAK'),  
5.           ('Walt','MIA'),  
6.           ('Buddy','ORD'),  
7.           ('Les','OMA')]  
8. #他们都要到美国的纽约集合，这是旅行的目的地  
9. #LGA是纽约的机场  
10. destination = 'LGA'  

从他们所在地到纽约是乘坐飞机，但是飞机有不同的航班，航班的价格也不一样。书中为我们准备了一份航班列表：schedule.txt，供我们使用，我节选其中一部分进行讲解,如下所示

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

他们用逗号分隔，依次是始发地、目的地、飞行时间、到达时间和价格。由于之前我们已经说了，这个旅游就是一天的时间的，一天之内往返，所以航班信息只有时间，没有日期。

我们必然会在代码中会使用这些航班信息，因为我们要针对每一个人，选出一对航班，让他从当前位置飞到纽约，再从纽约飞回他的家。所以，先让我们把这份文件读进代码，我们决定放在一个字典内，再次强调，字典就是键值对。
其中setdefault的函数，比较难理解，可以参考这里
代码如下：

[python] view plain copy

1. #我们将所有的航班信息读到一个字典内，以起止点为键，其他的为值  
2. flights={}  
3.   
4.   
5. for line in file('schedule.txt'):  
6.     origin,dest,depart,arrive,price = line.strip().split(',')  
7.   
8.   
9.     #其中setdefault是作为字典类的一个方法，主要的作用是应对一种情况：当一个键对应多个值,每一个值是一个元组，多个元祖组成一个列表，  
10.     #说白了，一个键对应一个列表，列表内有许多的元祖，一个元祖代表一个航班信息  
11.     #setdefault的含义是："如果没有，就设置"，如果有，就添加。  
12.     flights.setdefault((origin,dest),[])  
13.     #现在我们一行一行的把航班信息加入进去  
14.     #我认为，我们是对一个列表作为值，所以最外面是一个[]  
15.     #里面跟了一个元祖，这个元祖就是键  
16.     flights[(origin,dest)].append((depart,arrive,int(price)))  
17.     #通过print flight，我们可以看出在该字典中具体的存储情况。一个起止点作为一个值，对应了一个列表，而列表中有许多元祖，每一个元祖都是一个航班信息的起飞时间、降落时间和价格  

此时，我们需要加入一个函数，函数是将几点几分，转换为从凌晨0点开始经过了多少分钟，需要这个函数的原因：几点几分的比较非常不便，而单纯的分钟数，使得飞机的飞行时间、候机时间变得容易计算。

代码：

[python] view plain copy

1. def getminutes(t):  
2.     #strptime是将特定的时间格式转化为元组。  
3.     x=time.strptime(t,'%H:%M')  
4.     return x[3]*60+x[4]  
5. #我觉得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下表示。


显然，如果我们每次都这么读这一串数字是非常辛苦的，所以我们决定写一个函数，这个函数接受这一串数字列表，然后将人名、和起飞时间，起飞地点等信息都完整的打印出来。
代码：

[python] view plain copy

1. def printschedule(r):  
2.     for d in range(len(r)/2):#针对每一个人打印两个航班信息，我们只用重复r的一半次  
3.         name=people[d][0]  
4.         origin= people[d][1]  
5.         out=flights[(origin,destination)][r[2*d]]  
6.         ret=flights[(destination,origin)][r[2*d+1]]  
7.         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美元的罚款。

上述这是书中作者找到的关键的成本因素和金额度量。我们按照自己觉得合理的方式随意改动。

代码如下：

[python] view plain copy

1. def schedulecost(sol):  
2.     totalprice=0  
3.     latestarrival=0#最晚到底时间  
4.     earliestdep=24*60#最早离开时间，现在24*60是最好的情况，等一下会根据实际飞机情况发生改变  
5.   
6.   
7.     for d in range(len(sol)/2):  
8.         #得到每一个人的两次航班的价格并且加入到总价格中  
9.         origin = people[d][1]  
10.         outbound = flights[(origin,destination)][int(sol[2*d])]  
11.         returnf = flights[(destination,origin)][int(sol[2*d+1])]  
12.   
13.   
14.         #把钱加入到总价格里面去  
15.         totalprice+=outbound[2]  
16.         totalprice+=returnf[2]  
17.   
18.   
19.         #根据实际情况改变最晚到达时间和最早离开时间  
20.         if latestarrival<getminutes(outbound[1]):latestarrival=getminutes(outbound[1])  
21.         if earliestdep>getminutes(returnf[0]):earliestdep=getminutes(returnf[0])  
22.   
23.   
24.     #每一个人必须在机场等待直到最后一个来了才能出发  
25.     #每一个人必须在旅游结束时，为了最早离开的人能够赶上飞机，而来到机场等候  
26.     totalwait=0  
27.     for d in range(len(sol)/2):  
28.         origin = people[d][1]  
29.         outbound = flights[(origin,destination)][int(sol[2*d])]  
30.         returnf = flights[(destination,origin)][int(sol[2*d+1])]          
31.   
32.   
33.         totalwait+=latestarrival-getminutes(outbound[1])  
34.         totalwait+=earliestdep-getminutes(returnf[0])  
35.   
36.   
37.     if latestarrival<earliestdep: totalprice+=50  
38.     #如果这样看，在机场每多等一分钟就是1美元，而租车如果超过24个小时就罚款50美元  
39.     return totalprice+totalwait  

现在我成功打印出成本了，现在目标就是找出正确的数字序列，使成本最低。根据书中解释，一共是12次航班，每次有10种不同的航班。那么可以得到10的12次方个组合，就大约1000亿。如果我们将这1000亿个逐一比较，肯定能找到最佳答案。但是在计算机上耗费的时间很长。

随机搜索

所以我们需要求最优解。在正在开始学习最优化算法之前，让我们先来学一个叫随机搜索的东西。
随机搜索不是最优化算法，但是我们用它来评估其他算法的优劣。实际上，他就是随机产生结果。但是我们可以重复随机产生一堆结果...用来干什么我就不知道了。下面我们直接看代码吧，关键的部分我会在相应出解释。

代码：

[python] view plain copy

1. #domain代表随机产生的数字的个数和每一个数字的范围，是一个列表，列表里每个元素里面是一个元组，每个元组有2个元素，一个是上限，一个是下限  
2. #costf就是成本函数，  
3. #每次随机产生一组结果的时候，我们将会使用costf进行一下测试，看看效果如何  
4. def randomptimize(domain,costf):  
5.     best=999999999  
6.     bestr=None  
7.     for i in range(10000):#我们打算随机产生1000次结果，从这1000次结果中选择一个最好的  
8.         #很显然randint是产生在一定范围内的随机数,显然由于下一句右边等号里的for，将会产生一个循环  
9.         r=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]  
10.         cost=costf(r)  
11.   
12.   
13.         #每次得到成本我们都判断一次，如果更低，我们就置换  
14.         if cost<best:  
15.             best=cost  
16.             bestr=r  
17.     return bestr  
18.           

        
其中对domain这个理解需要比较到位，如下产生，打印出来看看就明白是什么了

[python] view plain copy

1. >>> domain=[(0,9)]*(len(people)*2)  
2. >>> print domain  
3. [(0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9), (0, 9),  
4.  (0, 9), (0, 9)]  

ok，在做了这么多准备工作之后，我们终于可以开始学习真在的求最优解的算法了。

爬山法

核心思维

通过随机搜索我们可以发现，每次寻找的解都是跳跃性的，前一个寻找的解和后一个寻找的解没有任何联系。所以在寻找最优解的过程中，效率不过，我们现在的想法是：让每一个解都对寻找最优解提供一点点帮助。爬山法的核心思维是：从一个随机解开始，然后在其临近的解集中寻找更好地题解，在本题中，就是更低的成本。核心思维有了，但是如何拿到临近解确实一个问题，或者说，如何设置临近解的产生，更加确切，在本题中，我们针对一个解，其中的每一个数字都加一或者减一，这样我们就能使得某一个人坐的飞机能够稍早一些或者稍晚一些。之所以加一或减一能够决定一个人稍稍早一些或者晚一些，这是因为：

经过查看flights{}字典发现，原来每一个键都是按照时间顺序排列的

如下所示，我记得我的代码并没有刻意为之排序，看样子是默认排序的。

[python] view plain copy

1. >>>   
2. {('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)]}  
3. >>>   

代码

核心思维结束了，让我们直接来看代码吧

[python] view plain copy

1. def hillclimb(domain,costf):  
2.     #先创建一个随机的解  
3.     sol=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]  
4.     while 1:#持续一个循环，直到在一次对每一个解减一或者加一之后没有任何改变时，就break  
5.         neighbors=[]  
6.         for j in range(len(domain)):  
7.             #解中的每一个元素都都会加一或者减一，加一产生一个解集，减一产生一个解集  
8.             if sol[j]>domain[j][0]:  
9.                 #如果很熟悉  
10.                 neighbors.append(sol[0:j]+[sol[j]-1]+sol[j+1:])  
11.             if sol[j]<domain[j][1]:  
12.                 neighbors.append(sol[0:j]+[sol[j]+1]+sol[j+1:])  
13.         #neighbors已经组装好了，现在在里面找最优解  
14.         current=costf(sol)  
15.         best=current  
16.         for j in range(len(neighbors)):  
17.             cost=costf(neighbors[j])  
18.             if cost<best:  
19.                 best=cost  
20.                 sol=neighbors[j]  
21.         if best == current:  
22.             break  
23.     return sol  

可以调用如下语句直接产生一个由爬山法制作出来的结果，当然请保证其他函数的存在。

[python] view plain copy

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

难点

其中有一个难点：sol[0:j]+[sol[j]-1]+sol[j+1:]

我做了个测试，看测试的还很清楚就能明白这个什么东西，加三个列表加起来，我觉得这个方法有点叼哇。最后一个错误请记住，要将三个列表相加的话，不要企图用数字。

[python] view plain copy

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

缺点

下面谈一下爬山法的缺点。我认为书本上的一幅图，非常好体现了爬山法缺点。简单来讲，爬山法只能找到局部最优解，不能找到全局最优解，而我们要找的当然就是全局最优解。当然我们可以重复使用爬山法，都使用不同的初始状态。

模拟退火算法

原理

受物理领域的一个启发。说白了很简单，根据刚刚所说的爬山法，我们很容易陷入局部最优解，这是因为每次我们都只允许成本小的存在（成为当前最优解，接着进行下一个循环）。我们没有给算法一个跳出这个局部的机会，所以，我们这个模拟退火算法的关键就是：如果成更高，我们让其有可能成为当前的最优解。关键的关键就是在可能的程度上，我们要求：刚开始的时候可能大，随着时间的推移（在模拟退火算法中是温度的降低，当温度降低到某个时刻），可能性越来越小，也就是只接受成本更低的数列。这就有了跳出局部最优解的机会。那么可能的程度，也就接受成本更大的概率，就使用一个公式来判断：

为什么刚开始的概率大呢？因为温度高，所以指数接近了0，所以概率几乎为1,随着温度的减少，高成本与低成本之间的差值将会越来越重要，差异越大，概率越低。这注定了：该算法只会接受稍稍高一点成本的解，而不会接受成本高出许许多多的解。

代码

具体执行代码如下：

[python] view plain copy

1. def annealingoptimize(domain,costf,T=10000.0,cool=0.98,step=1):  
2.     #和爬山法一样，先产生一个随机解，然后一切的改变都从这个随机解开始  
3.     vec=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]  
4.       
5.     while T>0.5:  
6.         #产生一个随机数，决定这次改变是改变数列中的哪一个随机数  
7.         i=random.randint(0,len(domain)-1)  
8.   
9.   
10.         #选择一个改变的方向，也就是说是增加还是减少  
11.         dir=random.randint(-step,step)  
12.   
13.   
14.         #复制随机解，然后对随机解进行改变，然后判断到底新的解好，还是后来产生的解好  
15.         vecb=vec[:]  
16.         vecb[i]+=dir  
17.         #这一段主要还是不让它超不过了最大最小值的限制  
18.         if vecb[i]<domain[i][0]:vecb[i]=domain[i][0]  
19.         elif vecb[i]>domain[i][1]:vecb[i]=domain[i][1]  
20.   
21.   
22.         #计算新产生的两次解的成本，然后对成本进行比较  
23.         ea=costf(vec)  
24.         eb=costf(vecb)  
25.           
26.         #or后面：表示接受更差的结果。仔细想想，原来概率的表示是如此完成的，注意前一个random()产生的数是在0到1之间。                    
27.         if(eb<ea or random.random()<pow(math.e,-(eb-ea)/T)):  
28.             vec=vecb  
29.   
30.   
31.         #没经过一次循环，改变温度，温度一改变，就会改变循环的次数和接受更差解的概率  
32.         #按一定比例降温  
33.         T=T*cool  
34.   
35.   
36.     return vec  

使用一下代码可以执行一次：

[python] view plain copy

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

经过我的验证，确实可以产生更不错的低成本，另外，偶尔也会得到一个较高成本的结果，记住一定要使用不同的参数多试一试。初始温度，冷却率，step的值等等。

遗传算法

原理

遗传算法是受生物学启发而产生了的。该类算法首先产生一个组解，称之为种群，对应种群里面的每一个解，我们都会计算其成本，显然一组解，应了一组成本。我们将成本排序。

如下图

依据这个排了序的结果，我们将会产生下一个子代。为了使下一个子代的成本更低，我们首先选出目前这一代的优良品种，由于排序，很容易选择前几名，我们可以也约定好，选择前多少名。

选了前几名之后，显然这一代的种群数量还不够，我们采用变异或者配对的方式来产生更多的子代，凑满种群的数量。

变异的做法就是指对一个解里面某一个小部分进行小的改变，如下图所示。6改为5,0改成1，当然这肯定是有实际含义的对吧。

配对又称交叉，我们将两个优良品种各取一部分，组成一个新的解。当然，实际上，要根据实际情况的不同而进行调整。如下图所示：

代码

然后我们凑够了下一代种群，之后又排序，重复上述过程，每一重复一次，我们为迭代一次，迭代的次数可以自己选定。
具体代码如下：

[python] view plain copy

1. #popsize：一个种群的规模大小  
2. #mutprob：种群中进行变异，而不进行配对的概率。  
3. #elite：在当前种群中被认为优秀的子代的比例，优秀的子代可以直接传入下一代  
4. #maxiter：迭代运行多少代  
5.   
6.   
7. def geneticoptimize(domain,costf,popsize=50,step=1,mutprob=0.2,elite=0.2,maxiter=100):  
8.     #方法中还在定义方法  
9.     #变异操作  
10.     def mutate(vec):  
11.         i=random.randint(0,len(domain)-1)  
12.         #完成第增加或减少的概率各一半  
13.         if random.random()<0.5 and vec[i]>domain[i][0]:  
14.             return vec[0:i]+[vec[i]-step]+vec[i+1:]  
15.         elif vec[i]<domain[i][1]:  
16.             return vec[0:i]+[vec[i]+step]+vec[i+1:]  
17.         else: return vec  
18.     #交叉操作：貌似用python编程是好快的说，我感觉比较复杂的句子只要两句么，还是我c/c++没学好  
19.     def crossover(r1,r2):  
20.         #为什么减2，其实想把这个一个数字列表划分为两段，再各取一半  
21.         i=random.randint(1,len(domain)-2)  
22.         return r1[0:i]+r2[i:]  
23.       
24. #构造初始种群  
25.     pop=[]  
26.     for i in range(popsize):  
27.         vec = [random.randint(domain[i][0],domain[i][1])for i in range(len(domain))]  
28.         pop.append(vec)  
29.           
30.     #每一代有多少优势物种，我们需要保留  
31.     topelite=int(elite*popsize)  
32.     #主循环  
33.     for i in range(maxiter):  
34.         #print pop #但是如果不加这句会使下一句出现一个bug，就是传过去的v是None,但是我讲pop全部打印出来的话，又没有问题  
35.         scores=[(costf(v),v) for v in pop]#列表里面，每一个元素都是一个元组，每一个元组是由一个数字和一个列表构成  
36.         scores.sort()  
37.         ranked=[v for (s,v) in scores]  
38.   
39.   
40.           
41.         #从中选择我们觉得优势的物种，然后保留  
42.         pop=ranked[0:topelite]  
43.   
44.   
45.         #如果种群数量不够，那么我们使用变异或者配对，产生新的后代个体  
46.         while len(pop)<popsize:  
47.             #变异的概率，这是由我们设定的,虽然这里是变异和配对只能选择其一，但是我认为是可以共同进行的  
48.             if random.random()<mutprob:#如果这样做，就是变异的少，交叉的多吧  
49.                 #变异  
50.                 c=random.randint(0,topelite)#注意是从优秀的子代中选出一个进行变异  
51.                 pop.append(mutate(ranked[c]))  
52.             else:  
53.                 c1=random.randint(0,topelite)#从优秀的子代中选择  
54.                 c2=random.randint(0,topelite)#从优秀的子代中选择  
55.                 pop.append(crossover(ranked[c1],ranked[c2]))  
56.               
57.         print scores[0][0]#注意打印的是成本          
58.                   
59.     return scores[0][1]#这里返回的是航班序列  

使用以下代码可以打印个结果:

[python] view plain copy

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

总结

到了这里算是一个阶段性结束。我们主要对一个模型抽象出了数字列表，然后对随机搜索、模拟退火算法、遗传算法进行了初步的理解和实现。我们发现在这里一系列问题中最困难的一步就是能否把问题潜在的解转化为数字列表。然后我个人觉得将整个问题转换为数字列表，或者能够运用编程的方式来解决就是很困难的。


此外，上述的求最优解的方法能发挥的功效呢？和问题本身有着密切的联系。上述求最优解的方法都依赖于一个事实：最优解应该接近去其他优解。然而，显示中可能存在这里的问题，如下图所示：


可以看出最优解的左右两边是否得陡峭，当产生了接近于最优解的解，我会认为其不是优解而排除，最后我们只能陷入了局部最优解之中，上图的左方。


如果放在航班例子中来看，就是我们如果从当天的第二次航班转到第三次航班时，比将其转到第八次航班更有可能降低成本。这是因为航班有序排列，对每一个解集的一个航班加一或者减一效果肯定比减5、加6来的效果好。但是如果航班处于无序状态，我们的求最优解的方式则不会有太大的作用，还不如随机搜索。其中的航班处于无序状态。什么意思？起飞的航班必须有序呀。


书中下面一节将了通过真实的航班搜索来完成，这个单独开一篇博客，因为情况涉及调用外部api。会稍稍复杂一些。

有偏好条件的求最优解

书中叫做涉及偏好的优化。
抽象性的描述就是：如果将有限的资源分配给多个表达了偏好的人，并尽可能使他们满意（或者根据他们的意愿，尽可能的满足他们的需要）。

这一次主要讲另一个例子，这个例子是学生宿舍分配求最优解问题

但是该代码可以很容易延伸到其他问题

• 在线纸牌游戏中的玩家牌桌分配
• 大型编程项目中开发人员的bug分配
• 家庭成员的家务分配

核心思想：从个体中提取信息，将其组合产生最优解。

学生宿舍例子

具体描述

学生宿舍例子的具体描述是：一共有5间宿舍，每间宿舍有两个隔间（显然美国人的宿舍是一人一个隔间），由10名学生来竞争住所。每个人学生有一个第一志愿和第二志愿。我们先用代码来描述一下基本情况，建一个py文件：dorm.py。代码如下：

[python] view plain copy

1. import math  
2. import random  
3. import MyOptimization  
4. #一个代表宿舍集合的列表，注意每一个宿舍有两个隔间  
5. dorms=['Zeus','Athena','Hercules','Bacchus','Pluto']  
6.   
7. #表述每个学生的第一志愿以及第二志愿  
8. prefs=[('Toby',('Bacchus','Hercules')),  
9.         ('Steve',('Zeus','Pluto')),  
10.         ('Andrea',('Athena','Zeus')),  
11.         ('Sarah',('Zeus','Pluto')),  
12.         ('Dave',('Athena','Bacchus')),          
13.         ('Jeff',('Hercules','Pluto')),  
14.         ('Fred',('Pluto','Athena')),  
15.         ('Suzie',('Bacchus','Hercules')),  
16.         ('Laura',('Bacchus','Hercules')),  
17.         ('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的范围确是：

[python] view plain copy

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

代码

下面就是打印出结果了，如果读懂了上述的思维，代码将会显得很容易理解了：

[python] view plain copy

1. #最关键的就是传入的vec这个序列的含义了  
2. #比如[0,0,0,0,0,0,0,0,0,0]  
3. #这并不代表某一个学生要选择dorms里面的第一个宿舍  
4. #它代表的含义是选择了剩余槽列表里面的第一个。  
5. def printsolution(vec):  
6.     slots=[]  
7.     #每个宿舍两个槽  
8.     for i in range(len(dorms)):slots+=[i,i]  
9.     print slots  
10.     #遍历每一个学生的安置情况  
11.     for i in range(len(vec)):  
12.         x = int(vec[i])  
13.         print slots  
14.         #从剩余槽中选择  
15.         dorm=dorms[slots[x]]  
16.         print prefs[i][0],dorm  
17.   
18.   
19.         #删除该槽  
20.         del slots[x]  

上述函数可以帮助我们打印出最终学生的选择，下面一个成本函数。成本函数的参数和printsolution是一I样的，我们针对每一个有效解（千万注意：该有效解不能直接看出某个选择选择了什么），计算了成本，如果成本低，我们就会保留。
成本函数：

[python] view plain copy

1. def dormcost(vec):  
2.     cost=0  
3.     #建立一个槽序列  
4.     slots=[]  
5.     #每个宿舍两个槽  
6.     for i in range(len(dorms)):slots+=[i,i]  
7.   
8.   
9.     #遍历每一个学生  
10.     for i in range(len(vec)):  
11.         x=int(vec[i])  
12.         dorm=dorms[slot[x]]  
13.         pref=prefs[i][1]  
14.         #首选成本值为0,次选成本值为1  
15.         if pref[0]==dorm:cost+=0  
16.         elif pref[1]==dorm:cost+=1  
17.         else cost+=3#不在选择之列则成本为3  
18.   
19.   
20.         #删除选中的槽  
21.         del slots[x]  
22.   
23.   
24.     return cost  

有了上面两个函数，我们只需要像之前解决航班问题一样，将domain和成本函数传入到求最优解的函数或者是随机搜索就可以了。


经过我实际运行，确实产生了最优解，非常棒。
在文档开头加入：import MyOptimization


代码：

[python] view plain copy

1. #[(0,9),(0,8),(0,7),(0,8).....(0,0)]  
2. domain=[(0,(len(dorms)*2)-i-1) for i in range(0,len(dorms)*2)]  
3. sol=MyOptimization.geneticoptimize(domain,dormcost)  
4. 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.

请注意，在这里我们使用求最优解算法的地方就在于求出题解（各个人的坐标），而每一个题解都有成本，成本就是交叉数，我们要请出成本最低的题解，

代码

接着来让我们看一下求成本的代码（也就是交叉线）：

[python] view plain copy

1. def crosscount(v):  
2.     #将数字序列转化为一个person:(x,y)的字典  
3.       
4.     loc=dict([(people[i],(v[i*2],v[i*2+1])) for i in range(0,len(people))])  
5.     total=0  
6.       
7.   
8.   
9.     #遍历每一对连线  
10.     for i in range(len(links)):  
11.         for j in range(i+1,len(links)):  
12.   
13.   
14.             #获得坐标位置  
15.             (x1,y1),(x2,y2)=loc[links[i][0]],loc[links[i][1]]  
16.             (x3,y3),(x4,y4)=loc[links[j][0]],loc[links[j][1]]  
17.               
18.             #下列几个算式涉及一个数学公式。不要太在意其复杂的表达方式，它只是遵循了数学公式  
19.             den=(y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)  
20.               
21.               
22.             #如果两线平行，则den==0  
23.             if den==0:continue  
24.   
25.   
26.             #否则，ua与ub就是两条交叉线的分数值  
27.             #如果要显示成小数点那样的形式，必须要用float，这个很重要。这是与书上不同的地方  
28.             ua=float(((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3)))/float(den)  
29.             ub=float(((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3)))/float(den)  
30.               
31.             #如果两条线的分数值介于0和1之间，则两线彼此交叉  
32.             if ua>0 and ua<1 and ub>0 and ub<1:  
33.                 total+=1  
34.     return total  

上述与书中不同的主要就是那个float，不加float的话，产生的数是四舍五入之后的数。

接着，我们可以用如下代码产生一个优秀的题解：
注意引入MyOptimization.py

[python] view plain copy

1. #因为我们决定将我们的坐标系建立在一个400*400的像素图中，下面写10-370是为了留点边缘，  
2. #所以我们可以可以知道每一个坐标的x或者y的范围都是10到370，而每一个人有x和y坐标，所以要乘以2  
3. domain=[(10,370)]*(len(people)*2)  
4. sol=MyOptimization.geneticoptimize(domain,crosscount)  
5. print sol  

我们现在已经求出最优解了，我们还想把它画出来。


画出来必须使用一个库文件，说白了就画图的呗。库的名字叫做Python Imaging Library。安装过程很简单，直接去下载，然后对应自己python版本下载一个，当然也要注意自己的平台。
然后就是代码，其实代码部分非常简单（不过让自己写还真的未必写的出来）

[python] view plain copy

1. from PIL import Image,ImageDraw  
2. def drawnetwork(sol):  
3.     #建立image对象  
4.     img=Image.new('RGB',(400,400),(255,255,255))  
5.     draw=ImageDraw.Draw(img)  
6.   
7.   
8.     #建立一个字典，准备打印位置，和上面计算成本一个意思  
9.     pos=dict([(people[i],(sol[i*2],sol[i*2+1])) for i in range(0,len(people))])  
10.   
11.   
12.     #划线，也就是说不会刻意去画点，全是线。  
13.     for (a,b) in links:  
14.         draw.line((pos[a],pos[b]),fill=(255,0,0))  
15.   
16.   
17.     #把人物名字绘制出来，看着人名就知道哪里是个点了吧  
18.     for n,p in pos.items():  
19.         draw.text(p,n,(0,0,0))  
20.     img.show()  

由于每次产生的题解都是不同的，所以每次画出来的图都不一样，但是现在可以看出来确实非常有趣了。下面贴一幅产生的图


但是这个图很丑（实际上很不错，我产生了很多次才产生出这个图的），关键是

• 两线夹角非常小
• 两点间的间距很近

为了画出好看的图，我们想要杜绝上述两种情况，以一个为例，如果我们不想要两个点靠的太近，那么我们可以里面的每对点进行遍历，如果像素靠的太近了，我们就加一点成本。在本题中，如果两点靠近的50个像素。就多加点成本，具体加了多少请看代码。

代码是在crosscount函数中先添加了一部分（在最后return之前）：

[python] view plain copy

1. def crosscount(v):  
2.     #将数字序列转化为一个person:(x,y)的字典  
3.       
4.     loc=dict([(people[i],(v[i*2],v[i*2+1])) for i in range(0,len(people))])  
5.     total=0  
6.       
7.   
8.   
9.     #遍历每一对连线  
10.     for i in range(len(links)):  
11.         for j in range(i+1,len(links)):  
12.   
13.   
14.             #获得坐标位置  
15.             (x1,y1),(x2,y2)=loc[links[i][0]],loc[links[i][1]]  
16.             (x3,y3),(x4,y4)=loc[links[j][0]],loc[links[j][1]]  
17.               
18.             #下列几个算式涉及一个数学公式。不要太在意其复杂的表达方式，它只是遵循了数学公式  
19.             den=(y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)  
20.               
21.               
22.             #如果两线平行，则den==0  
23.             if den==0:continue  
24.   
25.   
26.             #否则，ua与ub就是两条交叉线的分数值  
27.             #如果要显示成小数点那样的形式，必须要用float，这个很重要。这是与书上不同的地方  
28.             ua=float(((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3)))/float(den)  
29.             ub=float(((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3)))/float(den)  
30.               
31.             #如果两条线的分数值介于0和1之间，则两线彼此交叉  
32.             if ua>0 and ua<1 and ub>0 and ub<1:  
33.                 total+=2  
34.     for i in range(len(people)):  
35.         for j in range(i+1,len(people)):  
36.             #获得两个结点的位置  
37.             (x1,y1),(x2,y2)=loc[people[i]],loc[people[j]]  
38.   
39.   
40.             #计算两结点之间的距离  
41.             dist=math.sqrt(math.pow(x1-x2,2)+math.pow(y1-y2,2))  
42.             #对间距小于50个像素的结点进去判罚。  
43.             if dist<150:total+=(1.0-(dist/150.0))  
44.     return total  

然后我也改了改其中的数据，产生了两幅还不错的图。

总结

算法的应用场合也许有很多，但是看你是否去用心思考。

所有的难点在于：

• 确定题解的表示法
• 成本函数的确定

实际上，我们在处理三个类型的例子的时候，优化算法其实只写了一个。但是针对不同的问题，费了很多周折才完成了成本函数和题解的确定。

具体来说，什么样的问题能够使用求最优解的算法来解决呢？书中提出两点

• 问题本身有一个定义好的成本函数
• 相似的解会产生相似的结果

虽然不是所有此类问题都可以用求最优解的算法来解决，但是我们都可以试着研究一下。

对项目的思考

我的项目，哪里能用最优化算法？如果我能想到，我觉得，我很厉害。哈哈

不过暂时没想的很完整。

会不会当我把几个推荐引擎设计好之后，针对某一个用户，求出不同推荐引擎所占的比重，这个时候会不会用一下求最优解的算法呢？

那恐怕要先回答成本函数怎么定？

这个感觉得依据大量用户的操作数据吧。我的意思收集了大量一个用户操作数据之后，就可以做。因为比如引擎A产生的歌曲，用户点了删除，那么说明了....什么？

好像有点复杂。再想想吧。

全部源代码

一个例子一个，共三份

旅游例子

[python] view plain copy

1. # -*- coding: cp936 -*-  
2. import time  
3. import random  
4. import math  
5. #将要去旅行的人，第一个是名字，第二个是目前所在地   
6. people = [('Seymour','BOS'),  
7.           ('Franny','DAL'),  
8.           ('Zooey','CAK'),  
9.           ('Walt','MIA'),  
10.           ('Buddy','ORD'),  
11.           ('Les','OMA')]  
12. #他们都要到美国的纽约集合，这是旅行的目的地  
13. #LGA是纽约的机场  
14. destination = 'LGA'  
15.   
16. #我们将所有的航班信息读到一个字典内，以起止点为键，其他的为值  
17. flights={}  
18.   
19. for line in file('schedule.txt'):  
20.     origin,dest,depart,arrive,price = line.strip().split(',')  
21.   
22.     #其中setdefault是作为字典类的一个方法，主要的作用是应对一种情况：当一个键对应多个值,每一个值是一个元组，多个元组组成一个列表，  
23.     #说白了，一个键对应一个列表，列表内有许多的元组，一个元组代表一个航班信息  
24.     #setdefault的含义是："如果没有，就设置"，如果有，就添加。  
25.     flights.setdefault((origin,dest),[])  
26.     #现在我们一行一行的把航班信息加入进去  
27.     #我认为，我们是对一个列表作为值，所以最外面是一个[]  
28.     #里面跟了一个元组，这个元组就是键  
29.     flights[(origin,dest)].append((depart,arrive,int(price)))  
30.     #通过print flight，我们可以看出在该字典中具体的存储情况。一个起止点作为一个值，对应了一个列表，而列表中有许多元组，每一个元组都是一个航班信息的起飞时间、降落时间和价格  
31.       
32.   
33.   
34. def getminutes(t):  
35.     #strptime是将特定的时间格式转化为元组。  
36.     x=time.strptime(t,'%H:%M')  
37.     return x[3]*60+x[4]  
38. #我觉得x[0]x[1]x[2]是年月  
39.   
40. def printschedule(r):  
41.     for d in range(len(r)/2):#针对每一个人打印两个航班信息，我们只用重复r的一半次  
42.         name=people[d][0]  
43.         origin= people[d][1]  
44.         out=flights[(origin,destination)][r[2*d]]  
45.         ret=flights[(destination,origin)][r[2*d+1]]  
46.         print '%10s%10s %5s-%5s $%3s %5s-%5s $%3s' % (name,origin,out[0],out[1],out[2],ret[0],ret[1],ret[2])  
47.   
48. def schedulecost(sol):  
49.     totalprice=0  
50.     latestarrival=0#最晚到底时间  
51.     earliestdep=24*60#最早离开时间，现在24*60是最好的情况，等一下会根据实际飞机情况发生改变  
52.     for d in range(len(sol)/2):  
53.         #得到每一个人的两次航班的价格并且加入到总价格中  
54.         origin = people[d][1]  
55.         outbound = flights[(origin,destination)][int(sol[2*d])]  
56.         returnf = flights[(destination,origin)][int(sol[2*d+1])]  
57.   
58.         #把钱加入到总价格里面去  
59.         totalprice+=outbound[2]  
60.         totalprice+=returnf[2]  
61.   
62.         #根据实际情况改变最晚到达时间和最早离开时间  
63.         if latestarrival<getminutes(outbound[1]):latestarrival=getminutes(outbound[1])  
64.         if earliestdep>getminutes(returnf[0]):earliestdep=getminutes(returnf[0])  
65.   
66.     #每一个人必须在机场等待直到最后一个来了才能出发  
67.     #每一个人必须在旅游结束时，为了最早离开的人能够赶上飞机，而来到机场等候  
68.     totalwait=0  
69.     for d in range(len(sol)/2):  
70.         origin = people[d][1]  
71.         outbound = flights[(origin,destination)][int(sol[2*d])]  
72.         returnf = flights[(destination,origin)][int(sol[2*d+1])]          
73.   
74.         totalwait+=latestarrival-getminutes(outbound[1])  
75.         totalwait+=getminutes(returnf[0])-earliestdep  
76.   
77.     if latestarrival<earliestdep: totalprice+=50  
78.     #如果这样看，在机场每多等一分钟就是1美元，而租车如果超过24个小时就罚款50美元  
79.     return totalprice+totalwait  
80.   
81. #domain代表随机产生的数字的个数和每一个数字的范围，是一个列表，列表里每个元素里面是一个元组，每个元组有2个元素，一个是上限，一个是下限  
82. #costf就是成本函数，  
83. #每次随机产生一组结果的时候，我们将会使用costf进行一下测试，看看效果如何  
84. def randomptimize(domain,costf):  
85.     best=999999999  
86.     bestr=None  
87.     for i in range(10000):#我们打算随机产生1000次结果，从这1000次结果中选择一个最好的  
88.         #很显然randint是产生在一定范围内的随机数,显然由于下一句右边等号里的for，将会产生一个循环  
89.         r=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]  
90.         cost=costf(r)  
91.   
92.         #每次得到成本我们都判断一次，如果更低，我们就置换  
93.         if cost<best:  
94.             best=cost  
95.             bestr=r  
96.     return bestr  
97.   
98. def hillclimb(domain,costf):  
99.     #先创建一个随机的解  
100.     sol=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]  
101.     while 1:#持续一个循环，直到在一次对每一个解减一或者加一之后没有任何改变时，就break  
102.         neighbors=[]  
103.         for j in range(len(domain)):  
104.             #解中的每一个元素都都会加一或者减一，加一产生一个解集，减一产生一个解集  
105.             if sol[j]>domain[j][0]:  
106.                 #如果很熟悉  
107.                 neighbors.append(sol[0:j]+[sol[j]-1]+sol[j+1:])  
108.             if sol[j]<domain[j][1]:  
109.                 neighbors.append(sol[0:j]+[sol[j]+1]+sol[j+1:])  
110.         #neighbors已经组装好了，现在在里面找最优解  
111.         current=costf(sol)  
112.         best=current  
113.         for j in range(len(neighbors)):  
114.             cost=costf(neighbors[j])  
115.             if cost<best:  
116.                 best=cost  
117.                 sol=neighbors[j]  
118.         if best == current:  
119.             break  
120.     return sol  
121.   
122. def annealingoptimize(domain,costf,T=10000.0,cool=0.98,step=1):  
123.     #和爬山法一样，先产生一个随机解，然后一切的改变都从这个随机解开始  
124.     vec=[random.randint(domain[i][0],domain[i][1])for i in range (len(domain))]  
125.       
126.     while T>0.5:  
127.         #产生一个随机数，决定这次改变是改变数列中的哪一个随机数  
128.         i=random.randint(0,len(domain)-1)  
129.   
130.         #选择一个改变的方向，也就是说是增加还是减少  
131.         dir=random.randint(-step,step)  
132.   
133.         #复制随机解，然后对随机解进行改变，然后判断到底新的解好，还是后来产生的解好  
134.         vecb=vec[:]  
135.         vecb[i]+=dir  
136.         #这一段主要还是不让它超不过了最大最小值的限制  
137.         if vecb[i]<domain[i][0]:vecb[i]=domain[i][0]  
138.         elif vecb[i]>domain[i][1]:vecb[i]=domain[i][1]  
139.   
140.         #计算新产生的两次解的成本，然后对成本进行比较  
141.         ea=costf(vec)  
142.         eb=costf(vecb)  
143.           
144.         #or后面：表示接受更差的结果。仔细想想，原来概率的表示是如此完成的，注意前一个random()产生的数是在0到1之间。                    
145.         if(eb<ea or random.random()<pow(math.e,-(eb-ea)/T)):  
146.             vec=vecb  
147.   
148.         #没经过一次循环，改变温度，温度一改变，就会改变循环的次数和接受更差解的概率  
149.         #按一定比例降温  
150.         T=T*cool  
151.   
152.     return vec  
153.   
154. #popsize：一个种群的规模大小  
155. #mutprob：种群中进行变异，而不进行配对的概率。  
156. #elite：在当前种群中被认为优秀的子代的比例，优秀的子代可以直接传入下一代  
157. #maxiter：迭代运行多少代  
158.   
159. def geneticoptimize(domain,costf,popsize=50,step=1,mutprob=0.2,elite=0.2,maxiter=100):  
160.     #方法中还在定义方法  
161.     #变异操作  
162.     def mutate(vec):  
163.         i=random.randint(0,len(domain)-1)  
164.         #完成第增加或减少的概率各一半  
165.         if random.random()<0.5 and vec[i]>domain[i][0]:  
166.             return vec[0:i]+[vec[i]-step]+vec[i+1:]  
167.         elif vec[i]<domain[i][1]:  
168.             return vec[0:i]+[vec[i]+step]+vec[i+1:]  
169.         else: return vec  
170.     #交叉操作：貌似用python编程是好快的说，我感觉比较复杂的句子只要两句么，还是我c/c++没学好  
171.     def crossover(r1,r2):  
172.         #为什么减2，其实想把这个一个数字列表划分为两段，再各取一半  
173.         i=random.randint(1,len(domain)-2)  
174.         return r1[0:i]+r2[i:]  
175.       
176. #构造初始种群  
177.     pop=[]  
178.     for i in range(popsize):  
179.         vec = [random.randint(domain[i][0],domain[i][1])for i in range(len(domain))]  
180.         pop.append(vec)  
181.           
182.     #每一代有多少优势物种，我们需要保留  
183.     topelite=int(elite*popsize)  
184.     #主循环  
185.     for i in range(maxiter):  
186.         #print pop #但是如果不加这句会使下一句出现一个bug，就是传过去的v是None,但是我讲pop全部打印出来的话，又没有问题  
187.         scores=[(costf(v),v) for v in pop]#列表里面，每一个元素都是一个元组，每一个元组是由一个数字和一个列表构成  
188.         scores.sort()  
189.         ranked=[v for (s,v) in scores]  
190.   
191.           
192.         #从中选择我们觉得优势的物种，然后保留  
193.         pop=ranked[0:topelite]  
194.   
195.         #如果种群数量不够，那么我们使用变异或者配对，产生新的后代个体  
196.         while len(pop)<popsize:  
197.             #变异的概率，这是由我们设定的,虽然这里是变异和配对只能选择其一，但是我认为是可以共同进行的  
198.             if random.random()<mutprob:#如果这样做，就是变异的少，交叉的多吧  
199.                 #变异  
200.                 c=random.randint(0,topelite)#注意是从优秀的子代中选出一个进行变异  
201.                 pop.append(mutate(ranked[c]))  
202.             else:  
203.                 c1=random.randint(0,topelite)#从优秀的子代中选择  
204.                 c2=random.randint(0,topelite)#从优秀的子代中选择  
205.                 pop.append(crossover(ranked[c1],ranked[c2]))  
206.               
207.         print scores[0][0]#注意打印的是成本          
208.                   
209.     return scores[0][1]#这里返回的是航班序列  
210.           
211.   
212. domain=[(0,9)]*(len(people)*2)  
213. #sol=geneticoptimize(domain,schedulecost)  
214. #printschedule(sol)  
215. #print schedulecost(sol)  
216.                   

宿舍分配例子

[python] view plain copy

1. # -*- coding: cp936 -*-  
2. import math  
3. import random  
4. import MyOptimization  
5. #一个代表宿舍集合的列表，注意每一个宿舍有两个隔间  
6. dorms=['Zeus','Athena','Hercules','Bacchus','Pluto']  
7.   
8. #表述每个学生的第一志愿以及第二志愿  
9. prefs=[('Toby',('Bacchus','Hercules')),  
10.         ('Steve',('Zeus','Pluto')),  
11.         ('Andrea',('Athena','Zeus')),  
12.         ('Sarah',('Zeus','Pluto')),  
13.         ('Dave',('Athena','Bacchus')),          
14.         ('Jeff',('Hercules','Pluto')),  
15.         ('Fred',('Pluto','Athena')),  
16.         ('Suzie',('Bacchus','Hercules')),  
17.         ('Laura',('Bacchus','Hercules')),  
18.         ('Neil',('Hercules','Athena'))]  
19.   
20. #最关键的就是传入的vec这个序列的含义了  
21. #比如[0,0,0,0,0,0,0,0,0,0]  
22. #这并不代表某一个学生要选择dorms里面的第一个宿舍  
23. #它代表的含义是选择了剩余槽列表里面的第一个。  
24. def printsolution(vec):  
25.     slots=[]  
26.     #每个宿舍两个槽  
27.     for i in range(len(dorms)):slots+=[i,i]  
28.     #遍历每一个学生的安置情况  
29.     for i in range(len(vec)):  
30.         x = int(vec[i])  
31.         #从剩余槽中选择  
32.         dorm=dorms[slots[x]]  
33.         print prefs[i][0],dorm  
34.   
35.         #删除该槽  
36.         del slots[x]  
37.   
38.   
39.   
40. def dormcost(vec):  
41.     cost=0  
42.     #建立一个槽序列  
43.     slots=[]  
44.     #每个宿舍两个槽  
45.     for i in range(len(dorms)):slots+=[i,i]  
46.   
47.     #遍历每一个学生  
48.     for i in range(len(vec)):  
49.         x=int(vec[i])  
50.         dorm=dorms[slots[x]]  
51.         pref=prefs[i][1]  
52.         #首选成本值为0,次选成本值为1  
53.         if pref[0]==dorm:cost+=0  
54.         elif pref[1]==dorm:cost+=1  
55.         else: cost+=3#不在选择之列则成本为3  
56.         #删除选中的槽  
57.         del slots[x]  
58.     print cost  
59.     return cost  
60.   
61. #[(0,9),(0,8),(0,7),(0,8).....(0,0)]  
62. domain=[(0,(len(dorms)*2)-i-1) for i in range(0,len(dorms)*2)]  
63. sol=MyOptimization.geneticoptimize(domain,dormcost)  
64. printsolution(sol)  

社交网络绘图

[python] view plain copy

1. # -*- coding: cp936 -*-  
2. import math  
3. import MyOptimization  
4. from PIL import Image,ImageDraw  
5.   
6. people=['Charlie','Augustus','Veruca','Violet','Mike','Joe','Willy','Miranda']  
7. links=[('Augustus','Willy'),  
8.        ('Mike','Joe'),  
9.        ('Miranda','Mike'),  
10.        ('Violet','Augustus'),  
11.        ('Miranda','Willy'),  
12.        ('Charlie','Mike'),  
13.        ('Veruca','Joe'),  
14.        ('Miranda','Augustus'),  
15.        ('Willy','Augustus'),  
16.        ('Joe','Charlie'),  
17.        ('Veruca','Augustus'),  
18.        ('Miranda','Joe')]  
19.   
20. def crosscount(v):  
21.     #将数字序列转化为一个person:(x,y)的字典  
22.       
23.     loc=dict([(people[i],(v[i*2],v[i*2+1])) for i in range(0,len(people))])  
24.     total=0  
25.       
26.   
27.     #遍历每一对连线  
28.     for i in range(len(links)):  
29.         for j in range(i+1,len(links)):  
30.   
31.             #获得坐标位置  
32.             (x1,y1),(x2,y2)=loc[links[i][0]],loc[links[i][1]]  
33.             (x3,y3),(x4,y4)=loc[links[j][0]],loc[links[j][1]]  
34.               
35.             #下列几个算式涉及一个数学公式。不要太在意其复杂的表达方式，它只是遵循了数学公式  
36.             den=(y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)  
37.               
38.               
39.             #如果两线平行，则den==0  
40.             if den==0:continue  
41.   
42.             #否则，ua与ub就是两条交叉线的分数值  
43.             #如果要显示成小数点那样的形式，必须要用float，这个很重要。这是与书上不同的地方  
44.             ua=float(((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3)))/float(den)  
45.             ub=float(((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3)))/float(den)  
46.               
47.             #如果两条线的分数值介于0和1之间，则两线彼此交叉  
48.             if ua>0 and ua<1 and ub>0 and ub<1:  
49.                 total+=2  
50.     for i in range(len(people)):  
51.         for j in range(i+1,len(people)):  
52.             #获得两个结点的位置  
53.             (x1,y1),(x2,y2)=loc[people[i]],loc[people[j]]  
54.   
55.             #计算两结点之间的距离  
56.             dist=math.sqrt(math.pow(x1-x2,2)+math.pow(y1-y2,2))  
57.             #对间距小于50个像素的结点进去判罚。  
58.             if dist<150:total+=(1.0-(dist/150.0))  
59.     return total  
60.   
61. def drawnetwork(sol):  
62.     #建立image对象  
63.     img=Image.new('RGB',(400,400),(255,255,255))  
64.     draw=ImageDraw.Draw(img)  
65.   
66.     #建立一个字典，准备打印位置，和上面计算成本一个意思  
67.     pos=dict([(people[i],(sol[i*2],sol[i*2+1])) for i in range(0,len(people))])  
68.   
69.     #划线，也就是说不会刻意去画点，全是线。  
70.     for (a,b) in links:  
71.         draw.line((pos[a],pos[b]),fill=(255,0,0))  
72.   
73.     #把人物名字绘制出来，看着人名就知道哪里是个点了吧  
74.     for n,p in pos.items():  
75.         draw.text(p,n,(0,0,0))  
76.     img.show()  
77.     img.save('未改进的图.jpg','JPEG')  
78.           
79. #因为我们决定将我们的坐标系建立在一个400*400的像素图中，下面写10-370是为了留点边缘，  
80. #所以我们可以可以知道每一个坐标的x或者y的范围都是10到370，而每一个人有x和y坐标，所以要乘以2  
81. domain=[(10,370)]*(len(people)*2)  
82. sol=MyOptimization.geneticoptimize(domain,crosscount)  
83. print drawnetwork(sol)