问题描述
广义旅行商问题：给定一个加权图G=(V, E,W)，其中图中的边上有权值（所有权值为正）。给定一个源顶点s，和一个目标顶点点t，以及一个包含k个顶点的查询集合Q={vi1,vi2,…,vik}，其中k<=5。要求设计一个动态规划算法计算从s到t并且经过集合Q中所有顶点的最短路径。
要求：
1. 利用自己熟悉的编程语言实现算法，并且从以下URL中下载图数据进行测试
http://snap.stanford.edu/data/roadNet-CA.html其中边的权值可以随机生成一个正数即可。
测试时，可以随机生成顶点s，t，和集合Q （集合Q中的元素可以随机生成3至5个顶点），要求计算最短路径长度以及输出最短路径，同时统计算法的时间消耗和空间消耗。K值越大，成绩越高
介绍
这道题目是本人在大三《算法设计与分析》期末大作业的选题~
个人在做这个算法大作业的时候呢，哎，好多坑要爬 QAQ，所以呢想和大家分享一下，我在做这个题目的时候的流程&想法
存在较大瑕疵的就是“验证算法的正确性”，说服力还是有点牵强~
旅行商问题 & 广义旅行商问题
旅行商问题( TSP )：即为在赋权图 G 上找一条费用最小的 Hamilton 回路, 即一条能够遍历图中的一切顶点, 而且起点与终点重合的路
广义旅行商问题(GTSP)： 中, 顶点集 V 变为 p 个点群( cluster ) 的并集, V=V 1 ∪V 2 ∪…∪V p , 目标是要找到一条能够遍历p 个点群的费用最小的 Hamilton 回路
思考
这道题在一开始看来会觉得有点不像GTSP问题，但是仔细想了一下又好像是了。
举个例子Q1->Q2的最短距离，中间经过了许多的点，而这些点可以看做是Q1，Q2的点集，Q1，Q2就相当于是个点群（其中不包含查询集合Q中的点(除了Q1，Q2)，且不包含S,T起点终点）。至于点群中的点，有没有重复经过我们不关心，我们只需要确保这些点群我们只经历过1次，这便足矣。
那么要做的有2部分。一，寻找最短路径 二，动态规划找最优解
过程
1.读取文件
下载图文件发现一共有200w个顶点，这相当于2^21量级。若使用矩阵来进行图的存储，会发现就算是使用char（1byte）来进行存储，也有2^21 * 2^21 = 2^42 （相当于4000 G）的大小。这条路行不通。
一开始构建图就开始遇到难题了，这要如何解决呢？
观察该文件一共有5,533,213 共有550w条弧（这里我把它当做是一个有向图了，不过还是维持了图的连通性，而且往返代价不一也挺符合现实的，算是误打误撞）
那就只能使用邻接表进行图的存储。
定义struct结构
struct Node{
char weight; //0-128节省空间
int nodeID; //连通点下标
Node *next;
};//大小为12byte则存储图的空间大小为
5.55 * 10^6 * 12 byte = 6.66 * 10^7 byte ≈ 66.6 MB左右
伪代码：
void loadGraph(char *filePath) {
initNodeList();
while(!endOfFile()){
read(start,end);
createRandomDataForNode();
nodeID = end;
insertintoNodeList(start);
}
}//reading completed使用邻接链表这就代表着之后在查找的时候会比较的久，好在该图每个点的出度平均只有2.5。查找不会太费时。
2.Dijkstra算法寻找最短路径
INITIALIZE
S=∅
Q=G.V
while Q≠∅
u = EXTRACT-MIN(Q); //查找点集V – S 中找到最短路径
S = S ∪ {u}
for each vertex v ∈ G.Adj[u]
RELAX(u,v,w)
//松弛操作：更新最短路径，其中需要增加一些判断条件以排除Q –{end}中的点。
//一旦找到了start->end的最短路径就退出。集合S：从源结点s到该集合中每个结点之间的最短路径已经被找到。
集合Q：存放所有的结点集合
EXTRACTMIN(Q)：算法不断的从结点集Q（V - S）中选择最短路径估计值最小的结点u，
S = S ∪ {u}： 并把u加入到集合S中。
RELAX（u，v，w）：对所有从u点发的边进行“松弛”
If v.d > u.d + w(u,v)
v.d = u.d + w(u,v)
v.piror = u;验证算法的正确性：
查看图文件发现0->5之间只有3条路径，刚好可以用来测试算法是否会找到最短路径
取开始点为0，结束点为5，进行简单的测试。经过Dijkstra算法计算以后，得到最短路径： 0 -> 1 -> 6 -> 5, 最小开销为177
控制台输出的Graph为图~ 以链表的形式展示 
按照输出的权值画出图（只画出了主要的部分，其余无关的省略……） 
在3条路径之中，算法选择了上面那一条，简单的验证了编写的代码没有出错。
49+31+97=177
测试性能： 
随便设置了一个开始点0，结束点1000，经过70s才算出了结果。
如果设置得稍微远一点的点，那么……，发现程序运行了1个小时也没得到结果。
考虑到之后K=5的时候，相当于有7个点（5 + 2（S,T）），找一对点就要花费1小时，甚至更多的时间，那程序运行效率实在太低。
不解决这个问题，那么无法进行下去
现在分析一下Dijkstra算法：
Dijkstra算法的时间复杂度为O（n^2）其中extraMin操作的复杂度为O(n)，是最为核心的基本操作，所以必须从这里下手。
查阅网上资料可以发现，使用二叉堆可以大幅度的提高性能。
在为一组数据寻找最优数据的时候我们通常都需要遍历全部数据，然后找到最大（最小）值，一旦问题的规模一增大，就凸显了该算法的低劣。
换种思维，我们先把这些数据排序好，然后再一下子找到我们需要的数据。但是排序算法最优时间复杂度为O(nlogn)，感觉没有起到多大的优化效果。
但是在计算机应用中，我们实际上很少需要对一堆固定的数据进行查找，多数情况下是我们查找的数据集在不断的变化（增添，删除，改值）。
如果我们在每次数据发生变化的时候，花费较少的代价去维护这组数据，使其有序，而每次查找都只需要从有序数据集中直接获取数据，那么优化目的就达到了。
取得最优值只需要O(1)的时间复杂度，而在堆中进行插入，删除，修改操作的事实，也只需要O(logn)的时间复杂度.
使用C++ stl中的 set，就可以实现这个功能。
Set容器的内部结构为红黑树，其内部的排序算法时间复杂度为O(nlogn)，而排序只需要在一开始的时候排一次就足够了。
之后的都是进行数据的维护而已。所以总体来讲其时间复杂度为O(logn)
使用这种办法，原本时间复杂度为O(n^2)的Dijikstra，现在只需要O（nlogn），在现在这道题目中，规模为200w，其优化的效果要多得多。log2 22000000 ≈20 / 2000000 = 1/10^5
先验证改良后的Dijkstra算法的正确性： 
20+7+62+57=146正解。
接下来对比一下改良后与改良前的差距，为方便起见，这里把改良前后命名为1.0 和 2.0
在规模较小的时候，拿测试正确性的图作为例子。
因为容器set的构建的时间复杂度是O(nlogn)，在规模较小的时候无法体现出其强大之处。
在增大规模看看，如果设定找寻0-10000之间的最短路径，对比一下1.0和2.0的差别（这里为了测试所以固定了开始点与结束点） 
70s 和1000s之间的差距 1分钟和15分钟之间的差距。
这是在VS的debug模式下运行的，所以时间看起来比较慢。之后用release模式，有编译器自带的优化，速度会快很多。
有了2.0的Dijkstra算法以后，就可以生成我们的邻接矩阵了，如前面所理解的，利用Dijkstra算法求出
的最短路径，并为此建立一个邻接矩阵存储。则在main函数部分则需要调用k*(k-1)+2k约等于k^2次Dijkstra算法。这样再次说明优化好了这个算法，实验就成功了一半。那么建立这个邻接矩阵算法的时间复杂度为
N为点的数量，而k为查询集合的大小那么相比之下n远大于k，则建立邻接矩阵的时间复杂度还是为O(nlogn)。编写好代码之后，运行得到邻接矩阵
到这里，就真正的看出来差距了。现在是随机生成的起点，终点，查询集合。
1.0建立距离表花费了1362.36s，而2.0版本则仅仅需要7s。从这以后都使用2.0版本，k=3作为例子。
1.0的时间复杂度为0(n^2)，2.0则为O(nlogn)
把这个邻接矩阵命名为SQT 
3. 动态规划求最优解
利用动态规划的思想，先从顶向下去寻找“要求最优解需要什么？”
我们可以得到一个状态树： 
{S,{Q1,Q2,Q3}}表示是S经过集合{Q1,Q2,Q3}中的点到达T终点的最短路径
那么我们要得到最上面的最优解就需要“从底向上”去寻找。既然要用到动态规划的方法，就需要用到一个动态表来对数据进行存储。如果使用蛮力法，k有多少，则就有k层循环，每层循环次数-1，则蛮力法的时间复杂度为O(n!)
表的大小需要多少？ 
对于k种组合，用2进制进行表示。拿k=3作为例子Q={Q1，Q2，Q3}
Q1 = 001, Q2 = 010, Q3 = 100
{Q2,Q1} = 011, {Q3,Q2} = 110, {Q3,Q1} = 101, {Q3,Q2,Q1} = 111 
开始填表：
动态表，（记录 i 经过 j 中的集合 到达T的距离） 
定义string * binaryJ = new string[(int) pow(2,k)], 用于记录动态表列（j）的二进制数据。
更改：考虑到倘若用这么一个大数组来存储，会耗费太多的内存，那么这里用时间去换空间，每次都去计算一次j的二进制值，编写一个integerToBinary()函数，得到的结果存放在binaryJ数组当中 int *binaryJ = new int[k]; 该变量在每次循环都会被赋予新值
使用2种循环进行填表则其时间复杂度为O(n * 2^n) ，列为最外层循环，一列列的对表进行填充。表的大小也是
n * 2^n
从第一列开始：
if(j == 0){
table[i][j] = matrix[i][k+1];
}填入2720,1906,1454（矩阵中Q1，Q2，Q3到T的距离）
第二列，这里Q1不能经过自己，则这里需要增加一个判断条件：
If(one_snum == 1 || (i == 0 && ones_num == 1))
//当j中‘1’的个数为1时，或i=0且‘1’个数为1时
onesPosition = GetOnesPosition(j);
if(onesPosition == i) //若‘1’的位置与行i相同，则说明该位置不填
else //若位置不同，则到矩阵中去获取数据更新表
matrix[i][onesPosition] + table[onesPosition][0];
//Q2 -> Q1 -> T = Q2->Q1 + Q1->T
获取其位置（从后到前），这是为了能够与邻接矩阵表相互对应，这样才能从矩阵中取得数据填充到表中。
那么第二列中应该填
S->T = S->Q1 + Q1->T = 4380 + 2720=7100
Q2->T = Q2->Q1 + Q1->T = 4067 + 2720=6787
Q3->T = Q3->Q1 + Q1->T = 3591+2720=6311。
第三列：
S->T = S->Q2 + Q2->T = 737+1906=2643
Q1->T = Q1->Q2 + Q2->T = 4414+1906 = 6320
Q3->T = Q3->Q2 + Q2->T = 1217+1906 = 3123
第四列：
S->{Q2,Q1}->T = S->Q2 + Q2->Q1->T / S->Q1 + Q1->Q2->T即
这个空填的为min{（S->Q2 + Q2->Q1->T）,( S->Q1 + Q1->Q2->T)}
此时j：011。研究一下规律
S->Q2在邻接矩阵（matrix）中的坐标为[i(0)，2]，Q2->Q1->T则为[2，001]，第一个和第二个‘2’为010中‘1’的位置
S->Q1在邻接矩阵（matrix）中的坐标为[i(0)，1], Q1->Q2->T则为[1, 010]，同理此处的‘1’为001中‘1’的位置
那么此处应填下min{737+6787=7524，4380+6320=10700} = 7524
Q3->{Q2,Q1}->T = min{ Q3->Q2 + Q2->Q1->T + Q3->Q1 + Q1->Q2->T }
=min{ matrix[i][010(2)] + table[2][001(1)]，matrix[i][001(1)] + table[1][010(2)] }
=min{1217+6787=8004，3591+6320=9911} = 8004
剩余的同理填入。要到最后一列我们才能正确（确认）的推出这条关系式
① 我们要把J给分离成2部分，j中有多少个‘1’那就有多少种组合（ons_num = groups），我们要在这些组合当初去筛选最小的值填入表中
② 怎么分离呢？应该分出“只有一个‘1’的”和“剩余的”的2部分，对于前者（只有一个‘1’的），我们要记录其‘1’的位置和其十进制值为多少，后者则只需要记录其十进制值
分离的细节：
一共有多少个1就分离多少组。建立一个symbol变量，因为binaryJ不能更改。Symbol的作用接下来揭晓。举个例子，当前的binaryJ = 01011。
定义g1,g2数组用来存放分组结果，g1为“前者”，g2为“后者”
我们还需要2个函数getLowest1Pos()，getHighest1Pos()来获取最高位‘1’和最低位‘1’
定义 l1Pos，h1Pos 存放获取的位置
循环：
h1Pos = getHighestPos(binary) = 4;
第一次：l1Pos = getLowest1Pos(symbol) = 1;
g1: 00001 g2:01010
Symbol[l1Pos] = 0; ------> symbol = 01010
第二次：l1Pos = getLowest1Pos(symbol) = 2;
g1: 00010 g2:01001
Symbol[l1Pos] = 0; ------> symbol = 01000
第三次：l1Pos = getLowest1Pos(symbol) = 4;
g1: 01000 g2:00011
Symbol[l1Pos] = 0; ------> symbol = 00000
此时l1Pos = h1Pos 则退出循环
可以看出每次G1只有l1Pos位置才置为1，G2相当于binaryJ中l1Pos位置为0
③ 前者的十进制值作为matrix的j，i为matrix的i（当前循环的i值）。‘1’的位置作为table的i。后者的十进制则只作为table的j值
④ 点集中若包含了 i 中相同的元素，则该位置不填入，
if(binary[k - i] == 1) //i-1是因为字符从0开始最后一列：
一共有3个1，则可以分为3组
100 011：matrix[i][3] + table[3][011(3)] = 1541+8004=9545
010 101: matrix[i][2] + table[2][101(5)] = 737+7539 =8276
001 110: matrix[i][1] + table[1][110(6)] = 4380+6828=11208
则该空填下8276
我们可以得到递推式子：
至此动态表填完了，最终结果为
,这个值就是我们所求的最短路径开销
但是突然又想起了一个问题，怎么去获得路径？
想到了个比较耗费内存的办法，定义一个
在每次进行筛选最小值的时候记录table的下标
举个例子
Table[1][0] = matrix[1][4]
String[1][0] = “”
Table[3][1] = matrix[3][1] + table[1][0]
String[3][1] = 10
Table[2][5] = matrix[2][3] + table[3][1]
String[2][5] = 31
Table[0][7] = matrix[0][2] + table[2][5]
String[0][7] = 25
在回溯寻找最短路径的时候，定义string minPath
int i = 0; int j = 7(2^k-1);
while(string[i][j] != “”){
i = string[i][j][0] – ‘0’;
j = string[i][j][1] – ‘0’;
minPath += (char)(i + ‘0’);
}初始 i = 0，j = 7;
第一次循环： i=2,j=5; minPath = “2”
第二次循环: i=3,j=1 minPath = “23”
第三次循环: i=1,j=0 minPath = “231”
但是写完好像才发现有个bug，这个坐标只有个位数，如果是十位数，百位数，那就比较麻烦了。
然后突发奇想看到“25”，“31”，“10”这个不就是一个int型数据么，只要25/10,25%10就可以提取出2和5出来
如果是十位数，百位数的话，只要定个数（把10换成100 or 1000，根据k来确定这个数）
再换种想法，不用10，100，1000.而用来
替代它
这样的话理论上可以测到k = 31，用unsigned int 的话k=32，但是同时也要考虑到动态表的大小 
额…………
k = 23的话 4 * 23 * 2^23 = 736 MB × 2
k = 24 1536 MB × 2
k = 25 3200 MB × 2
k = 26 6656 MB × 2
电脑内存是8G，那么看看把以上推导的理论转为代码之后运行的结果！能否跑到k = 25
验证正确性：
代码测试结果与上述推导得到的结果一致，均为8276~
总结一下伪代码：
While( j ++< 2^k){
onesNum = getOnesNum(j); //获取j的二进制中‘1’的个数
integer2Binary(binary,j); // j 转为二进制，存储在bianryJ中
while( i++ < k +1){
if( j == 0) table[i][j] = matrix[i][k+1]; //填充第一列
else if(onesNum == 1 || ( i== 0 && onesNum ==1){
onesPosition = getOnePosition(binary); //获取1的位置
if(onesPosition == i) //若‘1’的位置与行i相同，说明该位置不放置
else {
table[i][j] = payMatrix[i][onesPosition] +
table[onesPosition][0];
tmpPath[i][j] = onesPosition * k + 0;
}
}
else if(onesNum > 1 && binaryJ[k - i] != 1){
group = onesNum;
do{
divideInto2Group(g1,g2);
onePos = getOnesPosition(g1);
n = binary2Integer(g2);
min(minPay, matrix[i][onePos] + table[onePos][n]);
}while(group--)
Table[i][j] = minPay;
}
} 至此完成所有代码的编写。
我们还需要验证一个东西：
我们的Q1->Q2,……这些路径中不能存在查询集合中的其他元素！以及S,T点
为此特地编写函数，利用之前回溯路径的函数showPath()中进行修改，在回溯路径的过程中对S,Q,T出现的次数进行统计
S,T只作为起点和终点，所以出现次数为K，而查询集合Q既是起点又是终点，故出现2k次
数据测试
实验环境：intel 酷睿i7 内存8g
先是进行压力测试：k=24时，在x64平台下勉强跑完。k =25的话8g内存运行不了。如前面所计算的需要6g的闲置内存，尽管设备内存8g，但是空余的内存并没有这么多 
Dijkstra算法：
原因可能是：
set容器在每次进行增、删操作的时候，更新最大（小）堆的时候只有在最坏的情况才会全部进行排序，而大部分的情况是达不到logn的。而根据拟合的曲线来看，实际的Dijkstra算法趋于logn
动态规划算法：
从曲线的趋势上看来，这有点迷….，O(n * 2^n),有点不一样，在查询集合k=17的时候，动态规划的耗时才开始逐渐呈“类似”指数形式的增长。
尝试分析一波（猜想）：
17*2^17*4byte=8912896 byte= 8 MB。
8MB这个数字正好和我的3级缓存总和相差不远。所以前面3-17的阶段一直有高速缓存的帮助，才能运行的如此的快，但是随着规模的增大，cache无法容下如此大的数据，只好利用到内存，从此之后才呈现出理想的增长趋势 