集合上的動態規劃

賈樹丙發表於2013-09-02
動態規劃

例題：最優配對問題

　　空間裡有n個點P0,P1,...,Pn-1，你的任務是把他們配成n/2對（n是偶數），使得每個點恰好在一個點對中。所有點對中兩點的距離之和應儘量小。n<=20,|xi|,|yi|,|zi|<=1000.

樣例輸入：

20
1 2 3
1 1 1
5 6 2
4 7 8
2 3 1
1 4 7
2 5 8
3 6 9
1 2 5
2 3 6
4 5 2
7 8 5
4 5 1
-1 2 3
-1 -9 -7
0 0 0
100 0 0
9 5 1
7 5 3
5 5 5

分析：

　　既然每個點都要配對，很容易把問題看成如下多階段決策過程：先確定P0和誰配對，然後是P1，接下來是P2，……，最後是Pn-1.按照前面的思路，設d(i)表示把前i個點兩兩配對的最小距離和，然後考慮第i個點的決策——它和誰配對呢？假設它和點j配對（j< i），那麼接下來的問題應是“把前i-1個點中除了j之外的其它點兩兩配對”，它顯然無法用任何一個d值來刻畫——我們的狀態定義無法體現出“除了一些點之外”這樣的限制。

　　當發現狀態無法轉移後，常見的方法是增加維度，即增加新的因素，更細緻地描述狀態。既然剛才提到了“除了某些元素之外”，不妨把它作為狀態的一部分，設d(i,S)表示把前i個點中，位於集合S中的元素兩兩配對的最小距離，則狀態轉移方程為：

　　d(i,S) = min{|PiPj|+ d(i-1,S-{i}-{j}) | j屬於S}

　　其中|PiPj表示Pi和Pj之間的距離。邊界是d(-1,S )= 0.

記憶化搜尋的程式碼如下：

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cmath>
 4 #define maxn 20
 5 #define INF (1<<30)
 6 int n;
 7 struct Point{
 8     double x, y, z;
 9 } p[maxn];
10 double f[(1<<maxn)];
11 inline double dist(int i, int j)
12 {
13     return sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x) + (p[i].y-p[j].y)*(p[i].y-p[j].y) + (p[i].z-p[j].z)*(p[i].z-p[j].z));
14 }
15 inline double min(double x, double y)
16 {
17     if(x < y) return x; 
18     return y;
19 }
20 double d(int S){
21     //記憶化搜尋快得要逆天了！！（因為奇數的都沒有算！！ ）
22     if(f[S] != -1) return f[S];
23     int i, j;
24     double& ans = f[S];
25     ans = (S == 0? 0 : INF);
26     for(i = n-1; i >= 0; --i) if(S & (1<<i)) break;
27     for(j = 0; j < i; ++j) if(S & (1<<j))
28       ans = min(ans, d(S^(1<<i)^(1<<j)) + dist(i, j));
29     return ans;
30 }
31 
32 int main(){
33     int i;
34     scanf("%d", &n);
35     if(n%2) {printf("No answer!\n"); return 0;}
36     for(i = 0; i < n; ++i) scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
37     for(i = 0; i < (1 << n); ++i) f[i] = -1;
38     printf("%.3lf\n", d((1<<n)-1));
39     /*
40     int cnt = 0;
41     for(i = 0; i < (1 << n); ++i) if(f[i]!=-1) ++cnt;
42     printf("%d\n", cnt);
43     */
44     return 0;
45 }

動態規劃的程式碼如下：

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cmath>
 4 #define maxn 20
 5 #define INF (1<<30)
 6 int n;
 7 struct Point{
 8     double x, y, z; 
 9 } p[maxn];
10 double f[maxn][(1<<maxn)];   
11 inline double dist(int i, int j)
12 {
13     return sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x) + (p[i].y-p[j].y)*(p[i].y-p[j].y) + (p[i].z-p[j].z)*(p[i].z-p[j].z));
14 }
15 inline double min(double x, double y)
16 {
17     if(x < y) return x; 
18     return y;
19 }
20 int main(){
21     int i,j;
22     scanf("%d", &n);
23     if(n%2) {printf("No answer!\n"); return 0;}
24     for(i = 0; i < n; ++i) scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);   
25     for(i = 0; i < n; ++i)
26         for(int S = 0; S < (1 << (i+1)); ++S){
27             double& ans = f[i][S];
28             ans = (S==0 ? 0 : INF);
29             if(S & (1 << i)) {
30                 for(j = 0; j < i; ++j) if(S & (1 << j))
31                     ans = min(ans, f[i-1][S^(1<<i)^(1<<j)] + dist(i, j));
32             }
33             else if(i!=0) f[i][S] = f[i-1][S];  
34         }    
35         printf("%.3lf\n", f[n-1][(1<<n)-1]);
36         return 0;
37 }

 

隱含的階段

　　上面的方程可以進一步簡化。事實上，階段 i 根本不用儲存，它已經隱含在S中了——S中最大的元素就是i。這樣，可直接用d(S)表示“把S中的元素兩兩配對的最小距離和”，則狀態轉移方程為：

　　　　d(S)=min{|PiPj|+d(S-{i}-{j}) | j屬於S，i=max(S)}

狀態有2n個，每個狀態有O(n)種轉移方式，總時間複雜度為O(n2n). 

　　值得一提的是，不少人一直在用這樣的狀態轉移方程：

　　　　d(S)=min{|PiPj|+d(S-{i}-{j}) | i,j屬於S}

　　它和剛才的方程很類似，唯一不同是：i和j都是需要列舉的。這樣做雖然也沒錯，但每個狀態的轉移次數高達O(n2)，總時間複雜度為O(n22n)，比剛才的方法慢。這個例子說明：即使用相同的狀態描述，減少決策也很重要。

　　那麼如何求出S中的最大元素呢？用一個迴圈判斷即可。當S取遍{0,1,2,...,n-1}的所有子集時，平均判斷次數僅為2。

1 for(int S=0; S<(1<<n); S++)
2 {
3     int i,j;
4     d[S] = INF;
5     for(i=0; i<n; i++)
6         if(S && (1<<i)) break;
7     for(j=i+1; j<n; j++)
8         if(S & (1<<j)) d[S] <?= d[S^(1<<i)^(1<<j)] + dist(i, j));
9 }

　　注意，在上述的程式中求出的i是S中的最小元素，而不是最大元素，但這並不影響答案。另外，j的列舉只需從 i+1開始——既然 i 是S中最小元素，則說明其他元素自然均比i大。最後需要說明的是S的列舉順序。不難發現：如果S‘是S的真子集，則一定有S'<S,因此S遞增的順序計算，需要用到某個d值時，它一定已經算出來了。

　　如果S'是S的真子集，則一定有S'<S。在用遞推法實現子集動態規劃時，該規則往往可以確定計算部分。

優化後程式碼如下：

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cmath>
 4 #define maxn 22
 5 #define maxs (1<<maxn)
 6 #define INF 1e10
 7 int n; 
 8 double f[maxs], x[maxn], y[maxn], z[maxn];
 9 inline double min(double x, double y)
10 {
11     return x<y?x:y;
12 }
13 inline double dist(int i, int j)
14 {
15     return sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) + (z[i]-z[j])*(z[i]-z[j]));
16 }    
17 int main(){
18     scanf("%d", &n);
19     for(int i  = 0; i < n; ++i) scanf("%lf%lf%lf", x+i, y+i, z+i);
20     int U = (1 << n);
21     for(int S = 0; S < U; ++S){
22         int i;
23         for(i = n-1; i >= 0; --i) if(S & (1<<i)) break;
24         f[S] = (S==0) ? 0 : INF;
25         for(int j = 0; j < i ; ++j) if(S & (1<<j)){
26           f[S] = min(f[S], f[S^(1<<i)^(1<<j)] + dist(i, j));
27         }    
28       }
29     printf("%.3lf\n", f[U-1]);
30     return 0;
31 }

 

相關文章