《算法導論》思考題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;
}