fc2ブログ



長方形探索: 最大長方形の面積 Largest Rectangle

2008.03.28 20:47  動的計画法

長方形探索

n × n のマス目のそれぞれに 1 または 0 が記してあり、その中から 1 だけから成る最大の長方形の面積を求めて下さい。

これは前に考えた正方形探索の応用で、今回は最大の長方形を探します。この問題も正方形探索で用いたアルゴリズムを応用して動的計画法で解くことができますが、正方形探索ほど単純な式では解決できません。左上角から右下角に向かって個々の要素を計算していく過程で、既に計算された左と上の要素を利用していきますが、長方形探索の場合、W[i][j] の値はそこから左上方向に向かってできる正方形の辺の長さではなく、そこから左上方向に向かってできる全ての長方形の情報を記録する必要があります。このアルゴリズムの詳しい解説が、プログラム・プロムナードに掲載されています。このアルゴリズムは、現在の要素を求めるために重複を避けながら左と上の要素をマージしたりと、プログラムがやや複雑になります。

ここでは、先に考えたヒストグラム内の最大長方形の面積を求めるアルゴリズムを応用したアルゴリズムを示します。このアルゴリズムは上記の正方形探索を応用したアルゴリズムよりもプログラムが簡単になるだけでなく、より効率の良いアルゴリズムです。アルゴリズムはいたって簡単で、まず各要素について上に向かって 1 が何個連続しているかを示すテーブル T を作ります。次に、T の各行をヒストグラムの入力と見なしヒストグラムの最大長方形の面積を求めるアルゴリズムを適用します。

最大長方形


001 #include<stdio.h>
002 #include<iostream>
003 #include<stack>
004 #include<algorithm>
005 #define MAX 105
006 #define BLOCK 0
007 
008 using namespace std;
009 
010 struct Rectangle{ int height;  int pos; };
011 
012 int getLargestRectangle( int size, int buffer[]){
013     stack<Rectangle> S;
014     int maxv = 0;
015     buffer[size] = 0;
016     
017     for ( int i = 0; i <= size; i++ ){
018         Rectangle rect;
019         rect.height = buffer[i];
020         rect.pos = i;
021         if ( S.empty() ){
022             S.push( rect );
023         } else {
024             if ( S.top().height < rect.height ){
025                 S.push( rect );
026             } else if ( S.top().height > rect.height ) {
027                 int target = i;
028                 while ( !S.empty() && S.top().height >= rect.height){
029                     Rectangle pre = S.top(); S.pop();
030                     int area = pre.height * (i - pre.pos);
031                     maxv = max( maxv, area);
032                     target = pre.pos;
033                 }
034                 rect.pos = target;
035                 S.push(rect);
036             }
037         }
038     }
039     return maxv;
040 }
041 
042 int size;
043 int buffer[MAX][MAX];
044 int T[MAX][MAX];
045 
046 int getLargestRectangle(){
047     for ( int i = 0; i < size; i++ ){
048         for ( int j = 0; j < size; j++ ){
049             T[i][j] = 0;
050         }
051     }
052     
053     for ( int j = 0; j < size; j++ ){
054         int sequence = 0;
055         for ( int i = 0; i < size; i++ ){
056             if ( buffer[i][j] == BLOCK ){
057                 sequence = T[i][j] = 0;
058             } else {
059                 T[i][j] = ++sequence;
060             }
061         }
062     }
063     
064     int maxv = 0;
065     for ( int i = 0; i < size; i++ ){
066         maxv = max( maxv, getLargestRectangle( size, T[i] ) );
067     }
068     
069     return maxv;
070 }
071 
072 int main(void){
073     while( cin >> size, size ){
074         for ( int i = 0; i < size; i++ ){
075             for ( int j = 0; j < size; j++ ){
076                 cin >> buffer[i][j];
077             }
078         }
079         cout << getLargestRectangle() << " ";
080     }
081     
082     return 0;
083 }

n 行についてO(n)の計算量なので、これはO(n2)のアルゴリズムとなります。



スポンサーサイト



テーマ : プログラミング - ジャンル : コンピュータ

| コメント(0) | トラックバック(0) | ↑ページトップ |




ヒストグラム中の最大の長方形の面積

2008.03.04 21:14  動的計画法

ヒストグラム中の最大の長方形の面積


データの離散型分布を表すヒストグラム(柱状グラフ)は、長方形を共通の基線上に並べた多角形として描画されます。これらの長方形は同じ幅を持ちますが、異なる高さのものを含みます。ヒストグラムを表す多角形に含まれる最大の長方形の面積を求めて下さい。入力は各データに対応する長方形の高さの列に対応した n 個の要素からなる1次元配列とします。(source: ACM/ICPC University of Ulm Local Contest)

以下のように、2重ループによって長方形の範囲を選び、そこにできる最大の長方形の面積を計算して更新していけば、全体での最大の面積を求めることができます。これはO(n2)のアルゴリズムです。



001 int getRectangleAreaBF( int size, int buffer[] ){
002 int maxv = 0;
003 for ( int i = 0; i < size; i++ ){
004 for ( int j = i; j < size; j++ ){
005 int minh = INT_MAX;
006 for ( int k = i; k <= j; k++ ){
007 minh = min( minh, buffer[k] );
008 }
009 maxv = max( maxv, minh*(j-i+1));
010 }
011 }
012 return maxv;
013 }



この問題は動的計画法を用いればO(n)のアルゴリズムで解くことができます。この解法は後で考える長方形探索などに応用することができるので重要です。

このアルゴリズムは部分問題の"未解決の"解をスタックに記憶していく、やや巧妙なテクニックを用います。スタックには、長方形の情報を記録します。各長方形は長方形の高さ height とその左端の位置 pos の情報を持たせます。まずスタックを空にし、ヒストグラムの各長方形 rect を左から右へ(i が 0 から n-1 まで)向かって順番に以下の処理を行います。


  1. スタックが空の場合:
      スタックに rect を追加する

  2. スタックの頂点にある長方形の高さが rect の高さより低い場合:
      スタックに rect を追加する

  3. スタックの頂点にある長方形の高さが rect の高さと等しい場合: 
      何もしない

  4. スタックの頂点にある長方形の高さが rect の高さより高い場合:

    • スタックが空でなく、スタックの頂点にある長方形の高さが rect の高さ以上である限り、スタックから長方形を取り出し、その面積を計算し最大値を更新する。長方形の横の長さは現在の位置 i と記録されている左端の位置 pos から計算できる。
    • スタックに rect を追加する。ただし、rect の左端の位置 pos は最後にスタックから取り出した長方形の pos の値とする




以下にアルゴリズムの流れを示します。位置 i におけるスタックには (i - 1) を右端としてできる未完成の長方形のリストが記録されます。図の 1 に示すように最初の5ステップでスタックには4つの未完成の長方形のリストが記録されます。

スタックの応用


図の 2 で6番目のデータである高さ 2 の長方形を追加します。ここから図の 3,4,5 で、2 よりも高い長方形をスタックから順番に取り出し面積を計算します。最後に取り出された高さ 3 の長方形の位置は 1 なので、この現在追加しようとしている高さ 2 の長方形の位置を 1 とします。この時点でスタックには高さ 1 と 2 の未完成の長方形が記録されています。

以下はPOJの解答プログラムです。



001 // pku 2559 (Largest Rectangle in a Histogram)
002 // O(n) algorithm
003 #include<stdio.h>
004 #include<iostream>
005 #include<stack>
006 #include<algorithm>
007 #define MAX 110000
008 typedef long long llong;
009
010 using namespace std;
011
012 struct Rectangle{ llong height; int pos; };
013
014 llong getRectangleArea( int size, llong buffer[]){
015 stack<Rectangle> S;
016 llong maxv = 0;
017 buffer[size] = 0;
018
019 for ( int i = 0; i <= size; i++ ){
020 Rectangle rect;
021 rect.height = buffer[i];
022 rect.pos = i;
023 if ( S.empty() ){
024 S.push( rect );
025 } else {
026 if ( S.top().height < rect.height ){
027 S.push( rect );
028 } else if ( S.top().height > rect.height ) {
029 int target = i;
030 while ( !S.empty() && S.top().height >= rect.height){
031 Rectangle pre = S.top(); S.pop();
032 llong area = pre.height * (i - pre.pos);
033 maxv = max( maxv, area);
034 target = pre.pos;
035 }
036 rect.pos = target;
037 S.push(rect);
038 }
039 }
040 }
041
042 return maxv;
043 }
044
045 main(){
046 int size;
047 llong buffer[MAX+1];
048
049 while(1){
050 scanf("%d", &size);
051 if ( size == 0 ) break;
052 for ( int i = 0; i < size; i++ ) scanf("%lld", &buffer[i]);
053 cout << getRectangleArea( size, buffer ) << endl;
054 }
055 }





テーマ : プログラミング - ジャンル : コンピュータ

| コメント(2) | トラックバック(0) | ↑ページトップ |




最大正方形の面積 正方形探索

2008.03.03 15:40  動的計画法

最大正方形


縦に n 行、横に n 列並べられた n × n のタイルがあります。いくつかのタイルには印(汚れ)が着いています。印のついていない(綺麗な)タイルだけから成るような最大の正方形の辺の長さを求めて下さい。各タイルの大きさは 1 × 1 メートルとします。

力任せのアルゴリズム

タイルの数が小さければこの問題は以下に示すループの入れ子で解くことができます。

  • 2重ループによって候補となる正方形の左上角(i, j)を決める

  • ループによって(i,j)から右下方向にできる正方形の幅 width を決める

  • 2重ループによって決められた範囲が正方形であるかチェックし、そうならば maxWidth を更新する



  • 正方形探索


    このアルゴリズムを実装すると以下のようになります。プログラムが5重ループになるので、O(n5) のアルゴリズムとなります。



    001 int size;
    002 char G[MAX][MAX];
    003
    004 int getLargestSquareBF(){
    005 int maxWidth = 0;
    006 for ( int i = 0; i < size; i++ ){
    007 for ( int j = 0; j < size; j++ ){
    008 for ( int width = 1; width <= size - max(i, j); width++ ){
    009 bool isSquare = true;
    010 for ( int a = i; a < i + width; a++ ){
    011 for ( int b = j; b < j + width; b++ ){
    012 if ( G[a][b] == BLOCK ) {
    013 isSquare = false;
    014 break;
    015 }
    016 }
    017 if ( !isSquare ) break;
    018 }
    019 if ( isSquare ) maxWidth = max( maxWidth, width );
    020 else break;
    021 }
    022 }
    023 }
    024 return maxWidth;
    025 }




    このアルゴリズムは n × n のタイル上の全ての正方形を見落とすことなく調べます。しかし明らかに無駄な処理が発生します。プログラムで break しているように印がついたタイルを見つけた時点でチェックする処理を終了させ次の候補の処理に移る、というようにアルゴリズムを改善したとしても、与えられたタイルのほとんどに印が着いていない場合には上記改善も意味がなく、O(n5)のアルゴリズムに変わりはありません。


    動的計画法による解法

    この問題は動的計画法を用いればO(n2)のアルゴリズムで解くことができます。

    データ
    W[i][j] が (i, j) から左上に向かってできる最大の正方形の辺の長さ(タイルの数)、というような配列 W を用意します。

    初期化
    上端と左端のタイル(i, j が 0 の場所にあるタイル)について
    T[i][j] に印が着いているならば、 W[i][j] = 0
    そうでなければ、W[i][j] = 1 とします。

    処理
    i, j > 0 について
    W[i][j] の値は、W[i-1][j-1], W[i-1][j], W[i][j-1] のうち一番小さい値に 1 を加えたものになります。例えば、下図の局面においては、左隣、左上隣、上隣の要素のうち最小の値は 2 です。これは、現在の位置を右下とした正方形の辺の長さは 2 + 1 より大きくすることはできないことを示しています。

    正方形探索


    計算が左→右、上→下の順に流れていくので、W[i][j] を計算するとき左隣、左上隣、上隣の要素の値はすでに計算済みなので、これらの解を利用することができます。このアルゴリズムを実装すると以下のようになります。プログラムが2重ループになるので、O(n2) のアルゴリズムとなります。



    001 int size;
    002 char G[MAX][MAX];
    003
    004 int getLargestSquare(){
    005 int W[MAX][MAX];
    006
    007 for ( int i = 0; i < size; i++ ) {
    008 W[0][i] = (G[0][i] == BLOCK) ? 0 : 1;
    009 W[i][0] = (G[i][0] == BLOCK) ? 0 : 1;
    010 }
    011
    012 int maxWidth = 0;
    013 for ( int i = 1; i < size; i++ ){
    014 for ( int j = 1; j < size; j++ ){
    015 if ( G[i][j] == BLOCK ) {
    016 W[i][j] = 0;
    017 } else {
    018 W[i][j] = min(W[i-1][j-1], min(W[i-1][j], W[i][j-1])) + 1;
    019 maxWidth = max( maxWidth, W[i][j] );
    020 }
    021 }
    022 }
    023
    024 return maxWidth;
    025 }


    テーマ : プログラミング - ジャンル : コンピュータ

| コメント(0) | トラックバック(0) | ↑ページトップ |




最長増加部分列 Longest Increasing Subsequence

2008.02.26 17:14  動的計画法

最長増加部分列

東西に並んだ高さの異なるn本の木があります。あるキコリがこれらの木を、西から東の方向へ木の高さが小さい順に並ぶように、いくつかの木を切り倒すことにしました。木の数を最大にするにはどの木を切れ(残せ)ばようでしょうか?

このような単純な(くだらない)問題でも、考えられる組み合わせを全て調べようとすると、それは各木を選ぶか選ばないかの組み合わせになるので、計算量はO(2n)となってしまいます。

この問題は、与えられた数列の最長増加部分列 Longest Increasing Subsequence (LIS) を求めることに帰着します。

最長増加部分列とは、与えられた数列 S

S = a1, a2 , , , an

の増加部分列 ( すべてのi, j (i < j)について ai < aj を満たす部分列 )
の中で長さが最大のものをいいます。

例えば、

S = 4 1 6 2 8 5 7 3

の増加部分列には

4 8
1 6 8
4 6 7
1 6 7
.
.

などがありますが、長さが最大のものは

1 2 5 7

です。

動的計画法によって与えられた数列の最長増加部分列をみつけることができます。
ここでは、n を数列の長さとすれば O(n2) の効率のアルゴリズムを紹介します。

与えられた数列を T とし、インデックスが 1 から始まるとします。また
L を、L[i] が「T[1]からT[i]までの要素でT[i]を最後に選んだときの LIS の長さ」を示す配列とします。

L[0] を 0 に初期化します。また T[0] は 0 (通常は非常に小さい値)に初期化します。
index 012345678
T041628573
L 0

以下 L[i] ( i が 1 から n まで)を計算します。

L[i] の計算:
j が 0 から i - 1 について、T[j] < T[i] を満たしかつ L[j] が最大である j を k とすると
L[i] = L[k] + 1 とします。

以下、計算の流れです。

index 012345678
T041628573
L 0 1

index 012345678
T041628573
L 0 1 1

index 012345678
T041628573
L 0 1 1 2

index 012345678
T041628573
L 0 1 1 2 2

index 012345678
T041628573
L 0 1 1 2 2 3

index 012345678
T041628573
L 0 1 1 2 2 3 3

index 012345678
T041628573
L 0 1 1 2 2 3 3 4

index 012345678
T041628573
L 0 1 1 2 2 3 3 4 3

L の要素で最大のものが、与えられた数列の LIS となります。

また、実際のLISを求めるためには、
P を P[i] が「T[1]からT[i]までの要素でT[i]を最後に選んだときのLISの1つ前の要素のインデックス」を示す配列とし、
各 L[i] を計算するときに、
P[i] = k
とすれば P の要素を辿ることで実際の LIS を求めることができます。
例えば、上記の例で P は以下のようになるので、
index 012345678
T041628573
L 0 1 1 2 2 3 3 4 3
P -1 0 0 1 2 2 4 6 4
p = 7 (LIS の最後のインデックス)
p = P[p] = P[7] = 6
p = P[p] = P[6] = 4
p = P[p] = P[4] = 2
よって、インデックスの列が 2 4 6 7 なので対応する LIS は
1 2 5 7
となります。

001 #include<iostream>
002 #define MAX 2000
003 using namespace std;
004 
005 int size, length;
006 int T[MAX+1], L[MAX+1], P[MAX+1];
007 
008 void printLIS(int i){
009     if ( P[i] == -1 ) return;
010     printLIS(P[i]);
011     cout << T[i] << endl;
012 }
013 
014 int LIS(){
015     T[0] = 0;
016     L[0] = 0;
017     P[0] = -1;
018     for ( int i = 1; i <= size; i++ ){
019         int k = 0; // max index
020         for ( int j = 0; j <= i - 1; j++ ){
021             if ( T[j] < T[i] && L[j] > L[k] ){
022                 k = j;
023             }
024         }
025         L[i] = L[k] + 1;
026         P[i] = k;
027     }
028     
029     // print result
030     int maxv = 0;
031     int maxi;
032     for ( int i = 1; i <= size; i++ ){
033         if ( maxv <= L[i] ){
034             maxv = L[i];
035             maxi = i;
036         }
037     }
038     
039     cout << maxv << endl;
040     cout << "-" << endl;
041     printLIS(maxi);
042 }
043 
044 main(){
045     cin >> size;
046     for ( int i = 1; i <= size; i++ ) cin >> T[i];
047     LIS();
048 }

テーマ : プログラミング - ジャンル : コンピュータ

| コメント(0) | トラックバック(0) | ↑ページトップ |




フィボナッチ数 Fibonacci Number

2008.02.26 16:06  動的計画法

Time-Space complexity など、動的計画法と再帰による手法の違いを表すよい例がイタリアの数学者フィボナッチによって定義されたフィボナッチ数(Fibonacci number)の生成です。n 番目のフィボナッチ数 Fn

Fn = Fn-1 + Fn-2

と定義されます。

これはフィボナッチがウサギの増殖をモデル化するために定義した数です。フィボナッチは、もし1年目に1つがいのウサギがいれば、n年目に生まれるウサギのつがいの数は、その1年前と2年前のウサギのつがいの和に等しいと推量しました。従ってフィボナッチ数列の最初の数項は F0 = 0, F1 = 1 とすれば
F2 = 1
F3 = 2
F4 = 3
F5 = 5
F6 = 8
F7 = 13
F8 = 21
F9 = 34
F10 = 55
.
.
となります。この数列は、自然界の現象やアプリケーションにも数多く出現します。

フィボナッチ数列は再帰的な式で定義できるので、それを以下のように簡単な再帰的なプログラムとして記述することができます。



001 int fibonacci( int i ){
002 if ( i == 0 ) return 0;
003 if ( i == 1 ) return 1;
004 return fibonacci(i - 2) + fibonacci(i - 1);
005 }



しかし、単純な再帰によるプログラムの計算量は指数関数時間になってしまいます。例えば、フィボナッチ数列の第5項を再帰で求めると以下のようになります。

再帰によるフィボナッチ数の生成


このように、最終的にはF(n)の結果がF(0)とF(1)の呼び出し回数の和になってしまいます。


以下のように、動的計画法を用いればF(n)をO(n)の多項式時間の計算量で求めることができます。


001 int F[MAX];
002
003 void makeFibonacci(){
004 F[0] = 0;
005 F[1] = 1;
006 for ( int i = 2; i < MAX; i++ ){
007 F[i] = F[i-2] + F[i-1];
008 }
009 }




動的計画法では全ての数列の値をメモリに記憶します。このプログラムはフィボナッチ数を小さい方から生成し記録していくので、F(i)を計算するときにはF(i-1)とF(i-2)がすでに計算されているのでそれを有効利用することができます。このように、小さい問題の解をメモリに記憶することにより、指数関数時間を多項式時間に改善し、計算時間を飛躍的に減少させることができます。

テーマ : プログラミング - ジャンル : コンピュータ

| コメント(0) | トラックバック(0) | ↑ページトップ |




ナップザック問題 Knapsack Problem

2007.02.07 13:47  動的計画法

knapsack.gif


大きさ w と価値 v を持った品物が N 個あり、これらを大きさ W のナップザックに入れたいとします。このとき、大きさの合計が W を超えず、価値の合計が最大になるような品物の組み合わせを求めたい。これがナップザック問題です。各品物を「選択する」か「選択しない」かの組み合わせなので、厳密には0-1ナップザック問題ともいいます。(各品物が複数個ある場合は、0123ナップザック問題と呼ばれます)

この問題を力任せで解こうとすれば、N 個の品物を「選択する」か「選択しないか」の全組み合わせを全て調べることになるので、計算効率は O(2N) となります。N が数十個でも、実用的な時間では計算できません。

品物の大きさ w、ナップザックの大きさ W がともに整数であれば、ナップザック問題は動的計画法により O (NW) の効率で厳密解を求めることができます。

C[i][w] が、大きさ w のナップザックに品物を 0 から i 個目まで考慮した時の価値の合計の最大値、となるような (N+1) × (W+1) の2次元配列 C を準備します。

C[i][w] の値は、

(1) C[i - 1][w - 品物 i の重さ] + 品物 i の価値、 または
(2) C[i - 1][w]

の大きい方となります。ここで、(1)はこの時点で品物 i を選択する、(2)はこの時点で品物 i を選択しない、という処理に対応します。ただし、(1)の場合は、品物 i の重さが w を超えない場合に限ります。

例を見てみましょう。

knapsackTables.gif

ナップザックの大きさ w が 0、または品物の数 i が 0 個の場合は、価値の合計が 0 なので、w = 0 と i = 0 の部分を 0 に初期化します。
大きさ w が 1 のナップザックに、品物 1 を入れてみる。
しかし、品物 1 の大きさは 2 で、容量を超えるので、無条件で品物 1 は選択しない(できない)。
大きさ w が 2 のナップザックに、品物 1 を入れてみる。
ここで、品物 1 の大きさは 2 なので選択可能。
(1) 選択した場合は 0 + 4 = 4 (斜めの矢印)、
(2) 選択しない場合は 0 (縦の矢印),
なので、C[i][w] = C[1][2] に 4 を記録する。
大きさ w が 3, 4, 5 についても同様。
大きさ w が 1 のナップザックに、品物 2 を入れてみる。
しかし、品物 2 の大きさは 2 で、容量を超えるので、無条件で選択しない。
大きさ w が 2 ~ 5 のナップザックに、品物 2 を入れてみる。
C[2][4] に注目してみましょう。ここでは
(1) C[1][2] + 品物 2 の価値 (= 4 + 5)または
(2) C[1][4] (= 4)`
のいづれか大きいほうを選びます。
大きさ w が 1 ~ 5 のナップザックに、品物 3 を入れてみる。
同様に C[i][w] を計算する。
品物 3 の重さが 1 なので、斜めの矢印の幅が 1 になる。
大きさ w が 1, 2 のナップザックに、品物 4 を入れてみる。
しかし、品物 4 の大きさは 3 で、容量を超えるので、選択できない。
容量を超えるとは、斜めの矢印が配列からはみ出してしまうことに相当する。。
大きさ w が 3, 4, 5 のナップザックに、品物 4 を入れてみる。
今度は、大きさが 3 なので、矢印の幅は 3 。
C[N][W] が価値の最大値。
ここから、矢印を伝っていくと、選ぶべき品物を記録することができる。
すなわち、斜めの矢印が品物の選択を表している。


#include<iostream>
#include<vector>
#define MAX 100
#define DIAGONAL 1
#define LEFT 0

using namespace std;
/**
* Target 624, 10130,
*
*/

struct Item{
    int value, weight;
};

int n;
Item items[MAX];
int C[MAX+1][MAX+1], G[MAX+1][MAX+1];
int weight;

void compute( int &maxValue, vector<int> &path ){
    for ( int w = 0; w <= weight; w++ ) {
        C[0][w] = 0;
        G[0][w] = DIAGONAL;
    }
    
    for ( int i = 1; i <= n; i++ ) C[i][0] = 0;
    
    for ( int i = 1; i <= n; i++ ){
        for ( int w = 1; w <= weight; w++ ){
            if ( items[i].weight <= w ){
                if ( items[i].value + C[i-1][w - items[i].weight] > C[i-1][w] ){
                    C[i][w] = items[i].value + C[i-1][w - items[i].weight];
                    G[i][w] = DIAGONAL;
                }else{
                    C[i][w] = C[i-1][w];
                    G[i][w] = LEFT;
                }
            }else {
                C[i][w] = C[i-1][w];
                G[i][w] = LEFT;
            }
        }
    }
    
    maxValue = C[n][weight];
    
    path.clear();
    int w = weight;
    
    for ( int i = n; i >=1; i-- ){
        if ( G[i][w] == DIAGONAL ){
            path.push_back( i );
            w -= items[i].weight;
        }
    }
    
    reverse( path.begin(), path.end() );
}

void input(){
    cin >> n >> weight;
    for ( int i = 1; i <= n; i++ ){
        cin >> items[i].value >> items[i].weight;
    }
}

main(){
    int maxValue;
    vector<int> path;
    input();
    compute( maxValue, path);
    
    cout << "max " << maxValue << endl;
    for ( int i = 0; i < path.size(); i++ ) cout << path[i] << endl;
    
    for ( int i = 1; i <= n; i++ ){
        for ( int j = 1; j <= weight; j++ ) {
            printf("%5d", C[i][j]);
        }
        cout << endl;
    }
    
}

| コメント(2) | トラックバック(0) | ↑ページトップ |




動的計画法 Dynamic Programming

2006.06.27 19:33  動的計画法

問題を解く有効な手法の1つが、最終的な解を求めるために、その問題の小さな部分問題の解を用いることです。一般的に小さな問題はより容易に分析・解決できるものです。この手法の代表的なものに、分割統治法(Divide and Conquer)と動的計画法 (Dynamic Programming: DP)があります。

分割統治法では問題を2つ(またはそれ以上)に分け、それらをそれぞれ解き、元の問題を解くためにそれらの分割して解かれた問題を統合して利用する、という処理を再帰的に行います。マージソートが分割統治法のよい例です。

一方動的計画法は、小さな部分問題を計算して得られた解をメモリに記録し、さらに大きい問題を解くためにそれら記録された結果を有効に使うことによって問題を解きます。

どちらも重要なプログラミング手法ですが、動的計画法は多くの問題の最適な解(最大値や最小値)を求めることができ、様々な問題に応用できる実用的なアルゴリズムです。




DP は、システムエンジニアやプログラマの方にとっては、必ず習得すべきプログラミングテクニックの1つだと思います。しかし、与えられた問題に DP が適用できるか否かの判断は難しいので、ある程度の経験が必要です。つまり、小さな問題の解を利用することによって、大きな問題の解が正しく求めることができるか否かの判断が難しいのです。

DP のスキルを伸ばすためには、DP が適用できる典型的な問題を解いてみるのが良いと思います。以下に応用問題をあげます(解法バレ注意してください)。











テーマ : プログラミング - ジャンル : コンピュータ

| コメント(0) | トラックバック(0) | ↑ページトップ |




連鎖行列積 Matrix Chain Multiplication

2006.06.26 01:36  動的計画法

l × m の行列 A と、m × n の行列 B が与えられた時、AB の積 C

C[i, j] = Σmk=1A[i, k] * B[k, j] ( 1 ≦ il, 1 ≦ jn)

を満たす l × n の行列になります。

matrixChain1.gif

この計算では、 l × m × n 回の掛け算が必要です。では複数の行列の積(連鎖行列積)について考えてみましょう。

(M1M2M3...Mn) のように、1 ≦ i ≦ n について、行列 Miri の行と ri+1 の列をもつ(下図参照) n 個の行列の掛け算を考えます。

matrixChain2.gif

この条件の下では、これらの行列の積は様々な順番で計算することができます。例えば、左から右へ向かって計算すると

(...((M1M2)M3)...Mn)

と表すことができ、右から左に向かって計算すると

(M1...(Mn-2(Mn-1Mn))...)

と表すことができます。さらに

(M1..(M2(M3M4)M5)...Mn)

などの、様々な順番が可能です。

ここで注目すべきことは、勿論掛け算の計算結果(連鎖行列積)は同じですが、行列の掛け算の順番によって、"掛け算の総回数"が異なってきます。例えば、上記の連鎖行列積に必要とする掛け算の総回数を、いくつかのケースについて考えてみましょう(2つの行列の掛け算には上記式の標準的な方法を用いることとします)。

Case 1. 左から右へ向かって: 86 回

matrixChainCase1.gif

Case 2. 右から左に向かって: 69 回

matrixChainCase2.gif

Case 3. 最適な順番: 36 回

matrixChainCase3.gif

連鎖行列積では、隣り合う2つの行列の掛け算にどのような方法を用いたとしても、「どのような順番で」掛け算を行うかが、総演算回数に大きく影響します。Matrix Chain Multiplication (連鎖行列積) は、掛け算の総回数を最小にする最適な連鎖行列の掛け算の順番を求める問題です。

連鎖行列積を計算する最適な順番を直接求めるために、全ての可能な順番を調べてみるアルゴリズムは現実的ではありません。行列の個数が大きくなれば、計算量が指数的に膨大になることは明らかです。以下に示すように、この問題には Dynamic Programming (DP) (動的計画法) が適用できます。

まず、(M1M2) を計算するための方法(順番)はただ 1 通りであり、r1 × r2 × r3 回の掛け算が必要です。同様に、(M2M3) を計算するための方法もただ 1 通りであり、r2 × r3 × r4 回の掛け算が必要です。つまり、(Mn-1Mn) を計算するための方法はただ 1 通りであり、rn-1 × rn × rn+1 回の掛け算が必要です。そして、これらの掛け算の回数をコストとして表に記録しておきます。

次に、(M1M2M3)、(M2M3M4)、...(Mn-2Mn-1Mn) を計算するための最適な方法を求めます。

例えば、(M1M2M3) を計算するための最適な方法を求めるためには、(M1(M2M3)) と ((M1M2)M3)のコストを計算し、小さい方のコストを (M1M2M3) のコストとして表に記録します。この時、

(M1(M2M3)) のコスト = (M1) のコスト + (M2M3) のコスト + r1 × r2 × r4
((M1M2)M3) のコスト = (M1M2) のコスト + (M3) のコスト + r1 × r3 × r4

となります。ここで注目すべきことは、これらの計算に用いた (M1M2) 及び (M2M3) は再計算する必要がなく、表から参照することができるということです。また、1 ≦ in について、(Mi) のコストは 0 であることに注意します。

一般的な形式では、(MiMi+1...Mi+j) を計算するための最適な方法は、 i + 1 ≦ ki + j における (MiMi+1...Mk-1)(Mk...Mi+j) の中から最小のコストを見つけることによって決定します。

例えば、、(M1M2M3M4M5) (つまり、i = 1、j = 4 の場合)のコストは、

k = 2 (M1)(M2M3M4M5) のコスト = (M1) のコスト + (M2M3M4M5) のコスト + r1×r2×r6
k = 3 (M1M2)(M3M4M5) のコスト = (M1M2) のコスト + (M3M4M5) のコスト + r1×r3×r6
k = 4 (M1M2M3)(M4M5) のコスト = (M1M2M3) のコスト + (M4M5) のコスト + r1×r4×r6
k = 5 (M1M2M3M4)(M5) のコスト = (M1M2M3M4) のコスト + (M5) のコスト + r1×r5×r6

のうちの最小値となります。


実装に向けて、さらに具体的に見てみましょう。

まず、cost[][] を cost[i][i + j]> が (MiMi+1...Mi+j) を計算するための最適なコストを記録する配列とします。

さらに、best[][] を best[i][i + j] が (MiMi+1...Mi+j) を (MiMi+1...Mk-1)(Mk...Mi+j) によって与えるような k を記録する配列とします。

以下に、プログラム例を示します。

#include<iostream>
#include<stdio.h>
#include<algorithm>

using namespace std;

void output( int, int, int, int);

#define MAX 100
#define INFTY (1 << 21)

int SIZE;
int R[MAX + 2];
int C[MAX+1][MAX+1]; // cost table
int B[MAX+1][MAX+1]; // best table

void compute(){
    for ( int i = 1; i <= SIZE; i++ )
    for ( int j = 1; j <= SIZE; j++ ) C[i][j] = INFTY;
    
    for ( int i = 1; i <= SIZE; i++ ) C[i][i] = 0;
    
    for ( int j = 1; j <= SIZE - 1; j++ ){
        for ( int i = 1; i <= SIZE - j; i++ ){
            for ( int k = i + 1; k <= i + j; k++ ){
                // C[i][i+j] = min( C[i][i+j], 
                //                  C[i][k-1] + C[k][i+j] + R[i]*R[k]*R[i+j+1]);
                int cost = C[i][k-1] + C[k][i+j] + R[i]*R[k]*R[i+j+1];
                if ( cost < C[i][i+j] ){
                    C[i][i+j] =  cost;
                    B[i][i+j] = k;
                }
                
            }
        }
    }
}

void order( int i, int j ){
    if ( i == j ) cout << "M" << i;
    else{
        cout << "(";
        order( i, B[i][j]-1 ); order(B[i][j], j);
        cout << ")";
    }
}

void input(){
    cin >> SIZE;
    for ( int i = 1; i <= SIZE + 1; i++ ) cin >> R[i];
}

void output(){
    cout << "minimum cost = " << C[1][SIZE] << endl;
    order(1, SIZE);
    cout << endl;
}

main(){
    input();
    compute();
    output();
}




Sample Input

6
40 20 30 10 20 20 30

Sample Output

minimum cost = 36000
((M1(M2M3))((M4M5)M6))


テーマ : プログラミング - ジャンル : コンピュータ

| コメント(1) | トラックバック(0) | ↑ページトップ |