双调欧几里得最短巡游路线问题

一个人能走地多快，往往和他积累的知识是正相关的。——动态规划

双调欧几里得最短巡游路线问题:

前言：

​ 本文将先用一个很容易想到的方法——利用标记的递归——来解决这个问题，稍后将不用标记来递归地解决这个问题，随后将证明问题的最优子结构，并稍微修改不用标记的递归来高效地解决问题，再根据递归式的特征决定解决子问题的顺序，用循环完成自底向上的动态规划方案，最后将问题的解——每个点连接的点——展示出来。

描述

对输入的设计：

​ 输入一个序列M，这个序列包含了一个点的x坐标和y坐标，这个序列已经根据X的大小关系进行了排序。这个序列的首个的下标是0，最后一个元素的下标是n。用结构体描述，就是：

    std::vector<struct POS{int x, y;}>map={……};

​ 注意：map[i].x并不一定就是i。

一些定义：

​ 函数d(i,j)是点M[i],M[j]的欧拉距离。代码实现中，虽然不会直接用i、j获得欧拉距离，但是下面的通过传递点的方案也行:

inline double eulerDistance(POS p, POS q){
    return sqrt(pow(p.x - q.x, 2) + pow(p.y - q.y, 2));
}

​ 本文一般使用“点n”来描述M[n]。
​ 递归的层：如:

​ $$f(x)=f(x-1)+1，第一层是f(x)，第二层是f(x-1)，第三层是f(x-2)……$$

一些说明 :

​ 由于双调巡回的翻译是bitonic tours,在后续的函数命名中，我将使用它或者它的缩写bnbt来命名一个函数。

​ 本文使用的递归表达式是这样的

\[f(x) = \begin{cases} x = 0, 1\\ x > 0, f(x - 1) + 1 \end{cases} \]

​ 和常见的数学递归表达式有所不同，但是这是为了表达一些复杂的递归式。

​ 程序中DBMAX是一个宏定义，#define DBMAX 9e300*9e300，实际上double类型无法存储这么大的数字，所以实际上这个表示无穷大，输出它的结果是"inf"。SIZE也是一个宏定义，它规定了动态规划的解的数量。
代码里面的goto标签全部没有作用——没有一处使用了goto语句，所以不需要思考什么地方出现了跳转。这些标签全部都是用来方便debug的——gdb可以在标签处打断点。不把语句写一行，目的是为了方便查看程序的执行情况。

一些容易发现的性质：

​ 题目中规定了，每个点都有不同的x坐标。
​ 待求的路径包含两条路径，起点和终点分别是点n和点0，分别叫这两条路径为Q1,Q2，路径用序列表示。比如\(Q1=<0, 2, 5>, Q2 = <0, 1, 3, 4, 5>\)，表示\(点0,2,5依次连接组成了Q1;点0,1,3,4,5依次连接组成了Q2\)。
​ 性质A：\(点p\in区间(0, n), p \in Q1或者p \in Q2\)

方案一，利用标记的递归穷举法：

​ 根据性质A，我们可以用一个数组来标记区间(0,n)里每一个点是哪条路径的。用数组bool where[] ，来标记开区间(0,n)上的点是在Q1上还是Q2上。如果where[i] = true, 那么点i在Q1上，否则在Q2上。当递归调用到n时，就计算两条路径的长度和，之后比较取最小值即可。map.size()刚好比n多一，所以这里用map.size() - 1来代替n。调用的接口是bitonicTours。容易得出，这个方法的时间复杂度是:

\[\Theta((n-1)*2^{n-1})=\Theta(n*2^n) \]

namespace recursion{
    double btnt(std::vector<POS>& map, std::vector<bool> where, int cur = 1){
        if(cur + 1 == map.size()){
            int u = 0, d = 0;
            double r = 0;//eulerDistance(map[1], map[0]);
            for(int i = 1; i < cur; i++){
                if(where[i] == true){
                    r += eulerDistance(map[d], map[i]);
                    d = i;
                }
                else{
                    r += eulerDistance(map[u], map[i]);
                    u = i;
                }
            }
            r += eulerDistance(map[u], map[cur]) + eulerDistance(map[d], map[cur]);
            std::cout << __func__ << ":" << r << std::endl;
            return r;
        }
        else{
            double m = btnt(map, where, cur + 1);
            where[cur] = true;
            m = std::min(btnt(map, where, cur + 1), m);
            return m;
        }
    }
    double bitonicTours(std::vector<POS>& map){
        std::vector<bool> where(map.size(), false);
        btnt(map, where);
    }
};

方案二，不用标记的穷举法：

​ 一般想不用标记，就需要根据问题特征来设计递归式。
​ 我们发现，Q1和Q2上面相邻的两个点（或者说Q1,Q2上的一条边的两个端点）之间分别是另一条路径的点。比如说，输入序列长度是9，那么Q1和Q2一定从点0开始，点8结束。如果路径Q1是序列<0,1,3,8>,那么Q2一定是<0,2,4,5,6,7,8>,Q1上相邻的两点3、8之间的4、5、6、7是Q2上的点。这个性质命名为性质B。这个很重要，请先记住这个。
​ 我们根据性质B优化方案一，来完成这个不用标记的穷举法。
​ 接下来设计递归式。递归式的参数应该有几个？
​ 我们知道，递归要用参数标记一个更小的子问题。比如，归并排序使用两个参数来标明待排序的范围，二分查找法用两个参数标明了待搜索的范围。所以要用若干参数标明子问题。

方案一的优化：

先来试试一个参数的。

​ f(i)表示从点n到i过点0回到n的双调巡回最短路。
​ 从点n出发，经过[0,n)间的某点，到达0，回到n，根据这个，写下递归式的一部分：

\[f(i)=min_{k∈[0,i)}(f(k)+……) \]

​ 发现有些写不下去了。要算出最短距离，对于f(i)这一层递归，我们要用上前面一层的递归内容。根据这个我们尝试写完它，当然感觉不论怎么写还是不正确。

\[f(i)=min_{k∈[0,i)}(f(k)+前面一层递归调用的参数i到当前递归调用的i之间的欧拉距离……) \]

​ 虽然感觉不正确，但是已经用上了别的递归调用层了。前面一层递归的参数，它属于(i,n]，凭借这个难以确定答案（如果你十分了解程序的堆栈结构，比如参数的位置、变量的位置、栈的增长方向，那么就可以做到了，当然，这个东西尽量别动）。
​ 由于需要使用当前递归调用层的前面一层或后面一层，不妨试试用两个参数。

两个参数的方案

​ 函数是f(i,j)。但是两个参数应该是什么意思？
​ 如果使用i,j表示点i、点j在同一条路径上，那么f(i,j)表示过边i,j的、过点0的最短路径。那么遇到的问题是，另一条路径上，如何把区间(i,j)之间的点连接到i之前，j之后的点。由此推断，这样子不太可行。
​ 那如果用i和j表明两条路径上的点，那么函数f(i,j)表明从i出发经过点0回到j的最短路径长。如此一来如果i在前面，那么根据性质B，我看可以连接(i, j)之间的点到j所在的路径上。举个例子，i=1，j=4时，那么2,3,4和4在同一条路径上。这个i<j且j不断减少的过程，暂时叫做跟随吧，比如刚刚这个例子里，叫做j跟随i。如此，通过不同的i就可以把每个点分配到两条路径上。
​ 为了保持j跟随i，那么就要使i小于j。注意，接下来一部分全部都是在条件i<j的情况进行的，请记住这点。

​ i<j时，这种跟随要到什么时候停下来？很容易想到，i+1=j时。
​ 好，由此写下递归式的一部分。

\[f(i,j) = \begin{cases} i+1<j, f(i,j-1)+d(j,j-1)\\ i+1=j,...... \end{cases} \]

​ 当然，也可以这样

\[f(i,j) = \begin{cases} i+1<j, \sum_{k\in(i, j)}^{f(i,k)+d(k,k+1)}+f(i,i+1)\\ i+1=j,...... \end{cases} \]

​ 由于两者都是在\(i+1<j\)的条件下把区间\((i,j]\)中的点连接起来，所以这两个等价，读者一定要自己试试，看看这两个是否等价。所以选择哪个似乎都无所谓。但是，根据前者确定动态规划填数组的顺序更加方便一点，因为它凸显了f(i,j)和f(i,j-1)的关系。如果使用后者，显示的是稍微复杂一点的关系：f(i,j)和区间\((i,j)\)之间的点的关系。为了让递归式看起来简单一点，也为了更加方便地实现动态规划，所以之后全部使用前者来写表达式。

然后i+1=j时该怎么办？

​ 有两种选择：
​ 在[0,i)之间选择一个点k，把i连接到k，也可以把j连接到k。试一试哪一个可行。

• i连接到k，那么j连接到区间(i,k)之间的点，那么递归式应该是
  \(i+1=j, \{min_{k∈[0,i)}^{f(k,j)+d(i,k)}\)

  出现一个问题，如果选择的k刚好是i-1，那么j不能链接到区间(i,k)之间的点，因为这里面没有一个点了。那么之后的递归调用，又需要得到之前的递归参数了。暂时放放这个方案。

• j连接到k，那么i连接到(i,k)间的点，那么应该用这个递归式。注意，为了保持\(f(i,j)\)的参数\(i<j\),这里就不是用f(i,k)。
  \(f(i,j) = i+1=j, \{min_{k∈[0,i)}^{ f(k,i)+d(j,k)}\)

​ 这里是j链接到k，(k, i)的点直接和i相连即可。
​ 好，上述问题解决。

写下稍微完整一点的递归式。

\[f(i,j) = \begin{cases} i+1<j, f(i,j-1)+d(j,j-1)\\ i+1=j, \min_{k∈[0,i)}^{ f(k,i)+d(j,k))} \end{cases} \]

这个递归式可以保持每次递归的参数i,j都保持i<j吗？

​ i+1<j时，一直调用表达式$$f(i,j-1)+d(j,j-1)$$,
​ i+1=j时，参数k<i，推断这可以做到。

这个递归式还没有结束的条件，接下来就想办法加入这个条件。

​ 当i=0时，那么就只用考虑j向0移动了。

\[f(i,j) = \begin{cases} i = 0\begin{cases} i+1<j, f(i,j-1)+d(j,j-1)\\ i+1=j, d(j,0) \end{cases}\\ i > 0\begin{cases} i+1<j, f(i,j-1)+d(j,j-1)\\ i+1=j, \min_{k∈[0,i)}^{ f(k,i)+d(j,k))} \end{cases} \end{cases} \]

​ 当然，在i=0时把0带入i，得到它这和下面的递归式等价

\[f(i,j) = \begin{cases} i = 0\begin{cases} 1 < j, f(0, j - 1) + d(j, j - 1)\\ j = 1, d(1,0) \end{cases}\\ i > 0\begin{cases} i+1<j, f(i,j-1)+d(j,j-1)\\ i+1=j, \min_{k∈[0,i)}^{ f(k,i)+d(j,k))} \end{cases} \end{cases} \]

​ 好，在i<j时的递归式就完成了。

接下来将递归式完善

​ 我们想知道的是从n出发的最短巡回路线长，而这个表达式只能求i<j时的答案，所以这个递归式不完整。
​ 回想方案一，它用一个数组标明了每个点分别是哪一条路径。那么n出发的两条边，分别连接到了第一个标记true和第一个标记false的点，那么求解f(n,n),就需要去寻找第一个和n连接的点，如此就出现了参数i<j的情况，于是就可以用之前的递归式了。

\[f(i,j) = \begin{cases} i = j = n\begin{cases} \min_{k∈[0,n)}^{ f(k,j)+d(j,k)} \end{cases}\\ i < j \begin{cases} i = 0\begin{cases} 1 < j, f(0, j - 1) + d(j, j - 1)\\ j = 1, d(1,0) \end{cases}\\ i > 0\begin{cases} i+1<j, f(i,j-1)+d(j,j-1)\\ i+1=j, \min_{k∈[0,i)}^{ f(k,i)+d(j,k))} \end{cases} \end{cases} \end{cases} \]

​ 注意，只要i<j=n，之后的递归调用，f的参数j都不可能是n。令n=3，然后带入f(2,3)试试？这一点在最后重构解时会用。

​ 不难看出，时间复杂度是\(\Theta(2^{n-1})=\Theta(2^n)\)

​ 接下来把递归的代码写下。由于\(i=n和j=n\)时，这个递归明显是递归套递归，所以设计了一个调用的接口bitonicTours，来求解f(n,n)，它将使得i不断减少到0，每次减少都调用了\(i<j\)的情况的递归。其中的std::cout来显示每个i,j的解。

namespace recursion{
    double _bitonicTours(std::vector<POS>& map, int i, int j){
        double d;
        if(i == 0){
            if(i + 1 == j)
                d = eulerDistance(map[i], map[j]);
            else{
                d = 0;
                for(int k = j - 1; k >= 0; k--)
                    d += eulerDistance(map[k], map[k + 1]);
            }
        }
        else{
            if(i + 1 == j){
                d = DBMAX;
                for(int k = i - 1; k >= 0; k--)
                    d = std::min(_bitonicTours(map, k, i) + eulerDistance(map[j], map[k]), d);
            }
            else{//i + 1 < j
                d = 0;
                for(int k = j - 1; k > i; k--)
                    d += eulerDistance(map[k], map[k + 1]);
                d = d + _bitonicTours(map, i, i + 1);
            }
        }
        std::cout << i << ',' << j << ':' << d << std::endl;
        return d、;
    }
    double bitonicTours(std::vector<POS>& map){
        double d = DBMAX;
        int i = map.size() - 1, j = map.size() - 1;
        for(i -= 1; i >= 0; i--)
            d = std::min(_bitonicTours(map, i, j) + eulerDistance(map[i], map[j]), d);
        return d;
    }
};

​ 随便来一组数据，

传入数据
{{0, 0}, {1, 2}, {2, 0}, {3, 1}}
cout的输出：
0,1:2.23607
1,2:4.23607
0,2:4.47214
2,3:6.47214
0,1:2.23607
1,2:4.23607
1,3:5.65028
0,3:5.88635

​ 从(0,1),(1,2),(1,2)的两个相同的解中，发现这个问题可能具有最优子结构和重叠子问题，下面开始证明这个问题具有最优子结构。
​ i+1=j时，如果f(i,j)是最优解，那么f(k,i)是最优解
​ 用反证法：p⇒q等价于¬q⇒¬p
​ 如果从[0,i)之间选择了一个不是k的点w，f(w,i)比f(k,i)更优。由于f(i,j)即f(i,i+1)过点0的路径的长度，与[i+1,n]之间路径长度的和是一个解，但是f(i,i+1)过点0的路径上的点，不和(i+1,n]相连，所以f(i,j)的过点0的路径的长度与[i+1,n]之间路径长度没有关联，所以选择k不是最优解。
​ 证毕。

​ i+1<j的情况，可以类似地证明。
​ 接下来就可以修改代码，这里用数组m[i][j]表示f(i,j)，使用自顶向下的方法解决问题。

namespace memorize{
    double m[SIZE][SIZE];
    double _bitonicTours(std::vector<POS>& map, int i, int j){
        if(m[i][j] < DBMAX)
            return m[i][j];
        double d;
        if(i == 0){
            if(i + 1 == j)
                d = eulerDistance(map[i], map[j]);
            else{
                for(int k = j - 1; k >= 0; k--)
                    d += eulerDistance(map[k], map[k + 1]);
            }
        }
        else{
            if(i + 1 == j){
                d = DBMAX;
                for(int k = i - 1; k >= 0; k--)
                    d = std::min(_bitonicTours(map, k, i) + eulerDistance(map[j], map[k]), d);
            }
            else{//i + 1 < j
                d = 0;
                for(int k = j - 1; k > i; k--)
                    d += eulerDistance(map[k], map[k + 1]);
                d = d + _bitonicTours(map, i, i + 1);
            }
        }
        std::cout << i << ',' << j << ':' << d << std::endl;
        m[i][j] = d;
        return d;
    }
    double bitonicTours(std::vector<POS>& map){
        double d = DBMAX;
        int i, j;
        for(i = 0; i < map.size(); i++)
            for(j = 0; j < map.size(); j++)
                m[i][j] = DBMAX;
        i = map.size() - 1, j = map.size() - 1;
        for(i -= 1; i >= 0; i--)
            d = std::min(_bitonicTours(map, i, j) + eulerDistance(map[i], map[j]), d);
        return d;
    }
}

之后使用自底向上的方法来解决问题。
如何确定填数组f[][]的顺序？
可以从递归式来确定。
不论i是不是0，只要i+1=j，那么就需要用f[k,i]来求f[i,j]；如果i+1<j，那么就需要用f[i,j-1]求f[i,j]
也可以画一个二维的数组，根据递归式来连接各点。如图

好，可以从上往下求解，当然还有一种求解顺序。在接下来的代码里，把每个i+1=j时做的选择记录了，这样方便重构解。印象中很少有要直接得到路径的题目。但是根据递归式我们可以发现i+1=j时有边j、k，i+1<j时，有边j、j-1。前者需要记录，后者根本不需要记录，直接根据j输出即可。在i+1=j时，可以把选择的点，给记录下来，记录到数组c吧，c的含义就是choice。 其中reconstructor的作用是根据c输出真正的路径，这个放在后面进行解释。
写下代码：

namespace dynamic{
    double m[SIZE][SIZE];
    int c[SIZE][SIZE];
    int n;
    double bitonicTours(std::vector<POS>& map){
        if(map.size() == 0)
            return 0;
        else if(map.size() == 1)
            return 0;
        else if(map.size() == 2)
            return 2 * eulerDistance(map[0], map[1]);
        int i, j, k;
        n = map.size() - 1;边
        for(i = 0; i < n; i++)
            for(j = 0; j <= n; j++)
                m[i][j] = DBMAX;
        m[0][1] = eulerDistance(map[0], map[1]);
dynamic:
        for(j = 2; j <= n; j++)
            m[0][j] = m[0][j - 1] + eulerDistance(map[j - 1], map[j]);
        for(i = 1; i < n; i++){
            double t;
            m[i][i + 1] = DBMAX;
            for(k = 0; k < i; k++){
                t = m[k][i] + eulerDistance(map[i + 1], map[k]);
                if(t < m[i][i + 1]){
                    m[i][i + 1] = t;
                    c[i][i + 1] = k;
                }
            }
            for(j = i + 2; j <= n; j++){
                m[i][j] = m[i][j - 1] + eulerDistance(map[j], map[j - 1]);
            }
        }
finish:
        double r = DBMAX;
        int s;
        for(k = n - 1, j = n; k >= 0; k--){
            if(m[k][j] + eulerDistance(map[k], map[j]) < r){
                r = m[k][j] + eulerDistance(map[k], map[j]);
                s = k;
            }
        }
        reconstructor(s);
        return r;
    }
};

重构解

​ 接下来就是根据c进行解的重构，这里重构的解的形式是"i->j"，意思是最优解在i、j之间有一条边。建议先把测试用例所产生的数组c抄一份，根据递归式自己手动重构一次解。前文在总递归式中，提及只要i<j=n，那么递归式将不会出现f(i,j)的参数j是n的情况。根据递归式，我们可以直接写下重构的代码。

namespace dynamic{
    void reconstructor(int i, int j = n){
        if(j == n)
            std::cout << i << "->" << n << std::endl;
        if(i == 0){
            if(j == 1)
                std::cout << 0 << "->" << 1 << std::endl;
            else{
                std::cout << j - 1 << "->" << j << std::endl;
                reconstructor(i, j - 1);
            }
        }
        else{
            if(i + 1 == j){
                std::cout << c[i][j] << "->" << j << std::endl;
                reconstructor(c[i][j], i);
            }
            else{
                std::cout << j << "->" << j - 1<< std::endl;
                reconstructor(i, j - 1);
            }
        }

    }
};

​ 对于之前的测试用例，输出是

2->3
1->3
0->2
0->1

多提供几组测试用例：

{{0, 0}, {1, 6}, {2, 3}, {5, 2}, {6, 5}, {7, 1}, {8, 4}}    ->    大约是25.58
{{1, 1}, {2, 3}, {3, 1}}    ->    大约是6.47
{{1, 1}, {2, 3}, {3, 1}, {4, 2}}    ->    大约是7.89

本文允许参考，但是禁止转载。