《算法导论》思考题15-1 双调欧几里得旅行商问题（动态规划）

　　欧几里得旅行商问题 是对平面上给定的n个点确定一条连接各点的最短闭合旅程的问题。图a给出了7个点问题的解，这个问题的一般形式是NP完全的,故其解需要多于多项式的时间。

　　J.K.Bentley建议通过只考虑双调旅程来简化问题，这种旅程即为从最左点开始，严格从左到最右点，再严格地从最右点回到最左点。图b显示了同样的7个点的问题的最短双调路线，在这种情况下，多项式的时间的算法是有可能的。

　　描述一个确定最优双调路线的O（n^2）时间的算法。可以假设任何两点的x坐标都不相同。

---

解法：

　　读了很多遍这个题，也看了网上好几篇关于这个问题的博客，有很多一部分是错误的却没有更正，于是自己实践整理了一份这个问题的解法。首先最大的疑惑在于理解什么叫做双调（Bitonic），读了很多解释的词条大概了解了双调的用处，我所理解的双调在这道题之中的思路就是将整个闭合的路线一个点展开，因为题中要求的是从左端向右端扫描，再从右端扫描回来。那么不妨将最左端的点作为出发点开始进行计算（从右端完全一致，不再赘述）。

　　我们需要一个辅助的矩阵a[iLen+1][iLen+1]，和所有动态规划问题一样矩阵的大小都是问题的规模+1，整个的计算过程是从左端到右端，也就是我们计算的最短路经从左端向右端生长的过程，辅助矩阵a[i][j]指的是p[i]到p[1]和p[1]到p[i-1]“共通”的最短路径，计算过程只利用矩阵的下三角部分，所以使前一个索引小与后一个索引。首先，将iLen个点存储到一个结构体/类数组之中，用来存放所有的点的坐标，记作数组p，p1,p2之间的距离记作dist（p1,p2）。

核心思想：

1）（前提）当我们计算第i个点或者将它并入最短路径的时候，前1.2.3...i-1个点已经形成了一条最短路径。

2）若要计算a[i][j]，新加入的点p[j]的选择分支有三种，也就是路线规划“生长”情况有且只有三种：

　　（一）当j-1〉i时，根据第（1）条我们需要做的就是在辅助矩阵中已经形成的子最短路径加上新增的边（按定义必然是先添加在更长的那半部分路线）

　　　　　就有 a[i][j]=a[i][j-1]+dist(j-1,j);

　　（二）当j-1=i时，就是一次总结前面的子最短路径生成更长的新子最短路径的时候。

　　　　　表示为 a[i][j]=min{(a[k][j-1]+dist(k,j)) , (a[k-1][j-1])+dist(k-1,j)....}////k为遍历值

 　  （三）当j=i 时，将整个图形封闭起来需要最后的两条边，实质上是（二）过程的一个分支

　　　　　即  a[i][j]=min{(a[k][j-1]+dist(k,j)+dist(j-1.j)).....}/////////k为遍历值

#include<iostream>
#include<vector>
#include<cstdlib>
#include<cmath>

#define MAX_VALUE 99999

class Point{
public:
	Point(double px=0,double py=0)
		:x(px),y(py){}////////Constructor & Default to be Zero
	double x;
	double y;
};

double dist(Point p1,Point p2){
	return sqrt((p2.x-p1.x)*(p2.x-p1.x)+(p2.y-p1.y)*(p2.y-p1.y));
}

int main(){
	const int iLen=7;
	Point p[iLen+1];//// Pointer List
	p[1]=Point(0,6);
	p[2]=Point(1,0);
	p[3]=Point(2,3);
	p[4]=Point(5,4);
	p[5]=Point(6,1);
	p[6]=Point(7,5);
	p[7]=Point(8,2);
	//////////////To minimize the time of Debugging  We get parameters initialize
	double a[iLen+1][iLen+1]={}; ////null-filled
	for(int j=3;j<=iLen;j++){
		for(int i=1;i<=j-2;i++){
			a[i][j]=a[i][j-1]+dist(p[j-1],p[j]);
		}
		a[j-1][j]=MAX_VALUE;
		
		for(int i=1;i<=j-2;i++){
			a[j-1][j]=(a[i][j-1]+dist(p[i],p[j]))<a[j-1][j]?(a[i][j-1]+dist(p[i],p[j])):a[j-1][j];
		}

	}
	double dMin=MAX_VALUE;
	for(int i=2;i<=iLen-1;i++){
		double tmp=a[i][iLen]+dist(p[i],p[iLen])+dist(p[iLen-1],p[iLen]);
		if(tmp<dMin)dMin=tmp;
	}
	///////////dMin records the answer
	std::cout<<dMin<<std::endl;
	return 0;
}