GLPKによるパズル解法 ぬりかべ

キューブ王

 最初の記事(111119)でのコードは速度の面で不満があった。そこで、今回書き換えた(130322)。説明も更新した。

 

[はじめに]

ニコリの「ぬりかべ」パズルの求解コードを、GLPKMathProg言語で記述してみた。ぬりかべパズル、分断禁止という制約がある。この制約、MathProgor整数計画法)での記述が難しいといわれている。しかし、スリリンのループ拒否コードと同様なやり方が可能で、簡単に攻略できる。

 前回の時点のコードでは、10*10問題を解くにはかなり時間が掛かっていたが、今回のコードではかなりの割合で1秒以内に解けるようになった。

 

[ルール]

 ニコリの記事に従い、ルールを示しておく。

 

 

以下のルールに従って盤面のマスをぬりつぶします。

(1)   数字が入っているマスは黒マスにはなりません。

(2)   数字は、その数字が含まれる、黒マスによって分断されたところ(シマとよぶ)のマスの数です。すべてのシマには数字が1つずつ入っていなければなりません。

(3)   すべての黒マスはタテヨコにひとつながりになっていなければなりません。

(4)   黒マスが2x2マス以上のカタマリになってはいけません。

以上が、ルール(制約)であるらしいが、実際には以下の制約もある。

(5) 2つのシマが辺を接してはいけません。(2)の1部ともいえるが、あえて明示。

さてこの内、ルール(3)は、「分断禁止」の制約と呼ばれている。この制約を、整数計画法で記述するのは難しい。

 なお、右の解答例をみれば、白マスが2*2の塊になった部分はない。しかし問題によっては2*2の塊が現れる。また、数字の白マスは必ず白いシマの端に現れている。しかし問題によっては、端でない場所に現れることもある。

 

[名称など]

以下はWikipediaの記事からの抜粋である。

------------------------------------------

「ぬりかべ」はニコリでつけられた名前である。海外においてもNurikabeの名称が使われているが、Cell Structure(部屋の構築),Islands in the Stream(島と潮流)などと呼ばれることもある。世界パズル選手権ではじめて出題されたとき(第7回)にはLay Bricks(煉瓦を敷く)というタイトルであった 

--------------------------------------------

当方、ファイル名などにはbricksを使うものと決めた。また個人的にwhite islands in the black currentという呼称が気に入ったので、それを使用する場合もある。なお、black currentは英辞郎によれば、「黒潮」である。さらに2*2の黒マスの塊は、「2*2の渦」と呼ぶことにした。

 

[問題の記述、データ構造]

 この問題、個別の白いシマを識別するほうがよさそうである。個別の白いシマ、1つ数字が入っているマスがある。このマス位置(アドレス相当)をシマに与えればいい。

 

 

左が問題の定義。真ん中が、個別のシマの数字位置の識別番号を与えたもの。右は、個別のシマ領域内の各マスに識別番号を与えたもの。識別番号が与えられてない部分が、黒潮である。

 次に、個別のシマのマス数を数えて、指定された数字と同じかどうか確かめる必要がある。

 

 

マス数を数えるには、個別のシマに木構造を導入すればいい(右図)。つまり、有向線(矢印)を置いていき、全てのマスが指定番号位置に到達できるようにする。そして、末端(葉ともいう)の方からマス数を数えていけばいい。右の図がその例。この場合、個別のシマが閉じてない線となり、さらに指定番号位置が端にあるので、数え上げは簡単。しかし、実際にはもう少し複雑な図になる場合もあるので、その点注意。

 

 

その具体例を示した。合流点での数え上げ方に注意。

黒潮部分、1つの白いシマ部分と同様に黒マス数を数えあげていける。その際、どこに数えあげていくか指定がある方がいい。その辺りを含め、実際のダンプ結果を示しておく(少し加工)

 

左のiniが問題定義配列。黒が保障される位置の1つに-1を設定している。この位置に黒マス数が数えあげられる。真ん中は、白いシマのマス数を数えあげたもの。右は、黒マスの数を数えあげたもの。

 

[解法コード]

 まずは、名称系や座標系を決める必要がある。

 

 

問題の盤面が2*2の場合として説明する。上下左右にマスを1つ追加して4*4とする。個別のマスが白かどうかをw[,]とする。[1,1]の位置のマスだとw[1,1]である。一般にはw[r,c]とする。なお、余分に追加したマスの色はコウモリ色である。灰色ではない。黒でありさらに白であるという特別な色である。

個別の白いマスに数字を与えるのだが、それはwban[r,c]とする。

また、有向線(矢印)wr[r,c],wu[r,c],wr[r,c],wd[r,c]とする。個別のマスから出て行く矢印である。従来、縦の辺や横の辺にも座標系を定めていて、その辺を横切る矢印という考えで矢印に座標を与えていた。今回は、マスから出て行く矢印という考えで座標系を定めている。つまり、あるセル[r,c]から出て行く矢印はwu[r,c]などとしてある。その点注意。

解法コードの基本の部分を順次に示していく。

 

          param CB:=1;                                                                                 #

## +++++++++++ white islands in the Balck Current ++++++++++++++++

          param Nr integer;                                                             # number of rows

          param Nc integer;                                                             # number of columns

          param N:=Nr*Nc;

          param ini{1..Nr,1..Nc} integer;                                              # board to define problem

          param SN:=0.001;

          var w{0..Nr+1,0..Nc+1} binary;                                              # 1/0: white/not

          var b{0..Nr+1,0..Nc+1} binary;                                              # 1/0: black/not

## ------------ initialization ---------------------------

          var wban{0..Nr+1,0..Nc+1} >=0;                                            # work board for counting white cells

s.t. INI0{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}: wban[r,c]=ini[r,c];#

s.t. INI2{r in 0..Nr+1,c in 0..Nc+1:                                                  # outers cell: black and white

                         r=0||r=Nr+1||c=0||c=Nc+1}: b[r,c]+w[r,c]=2;

s.t. BLACK0{r in 1..Nr,c in 1..Nc}: b[r,c]=1-w[r,c];                            # black is "not white" on the board.

## ------------ inhibit vortexes caused by black current ---------------

s.t. VORT{r in 2..Nr,c in 2..Nc}: b[r-1,c-1]+b[r-1,c]+b[r,c-1]+b[r,c]<=3;                           # bB

 

 param CBは、C言語の”#define CB”に相当する。これが1だと、黒マス数を数えあげるための制約が活性化される。

さて、param Nrは問題の盤面の行数、param Ncは列数である。Param Nはマス数である。シマに与える識別番号の最大値でもある。

 Param ini[,]は、問題定義配列。[ルール]節の最初の図の左端の図のような問題を与える。具体的に与えた例はこのページ最後尾の問題を参照。

 Param SNは、0に近い正数。300倍しても1以下である。そういう値。

 変数w[,]は、盤面の個別マスの白/(1/0)を定めるもの。変数b[,]は、盤面の個別マスの黒/(1/0)を定めるもの。

 変数wban[,]は、個別のシマ毎に白いマス数を数えあげるためのもの。binaryではなく整数でもなく、実数である点に注意。整数なら wban[r,c]=n を表現する代わりに3次元のbbb[r,c,n]を導入して bbb[r,c,n]=1 とするやり方もある。しかし、ここでは2次元実数のwban[r,c]である。

 制約INI0は、数字が与えられていたマスに、その数字を識別番号として与えるもの。

 制約INI1は、数字が与えられていたマスは白マスであるというもの。[ルール]節のルール(1)に対応。

 制約INI2は、拡大したマス部分はコウモリ色(黒かつ白)マスであるというもの。ルールからではなく、コードを記述する都合のもの。

 制約BLACK0は、本来の盤面上では、白でなければ黒であるというもの。

 制約VORTは、黒マス(黒潮)2*2の渦を作ってはならないというもの。[ルール]節のルール(4)に対応。

 

## ------------ id numbers for white cells ------------------------------------

              var wid{0..Nr+1,0..Nc+1} >=0;                  # id number for white cells

s.t. WID00{r in 0..Nr+1,c in 0..Nc+1: r=0||r=Nr+1||c=0||c=Nc+1}: wid[r,c]=0;             # outer cells have not id number

s.t. WID01{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}:             wid[r,c]=Nc*(r-1)+c;                             # center cells have id number

s.t. WID02{r in 1..Nr,c in 1..Nc}: wid[r,c]<=N*w[r,c];                                            # black cells have not id number

 

s.t. WID038{r in 1..Nr,c in 2..Nc}: wid[r,c-1]<=N*(b[r,c]+b[r,c-1])+wid[r,c];                            # wW-->wid=wid

s.t. WID039{r in 1..Nr,c in 2..Nc}: wid[r,c  ]<=N*(b[r,c]+b[r,c-1])+wid[r,c-1];          

s.t. WID040{r in 2..Nr,c in 1..Nc}: wid[r-1,c]<=N*(b[r,c]+b[r-1,c])+wid[r,c];                            # w --> wid=wid

s.t. WID041{r in 2..Nr,c in 1..Nc}: wid[r,  c]<=N*(b[r,c]+b[r-1,c])+wid[r-1,c];                         # W

 

 変数wid[,]は、個別の白マスに識別番号を与えるためのもの。

 制約WID00は、拡大した部分のマスのwid[r,c]は0というもの。

 制約WID01は、数字指定のマス位置のwid[r,c]に、識別番号 Nc*(r-1)+c を与えるもの。

 制約WID02は、黒マス位置のwid[r,c]に0を与えるもの。この行がないと、所々の黒マスに変な番号が付く。

 以降の4つの制約、WID038,039,040,041は、2つの白マスが辺を接しているならば、2つの白マスのid番号は同じ、というもの。最初の2つは、横位置の2つ。最後の2つは、縦位置の2つ。[ルール]節のルール(5)に対応。

 

## ------------ arrows for counting white cells ---------------------------

              var wl{0..Nr+1,0..Nc+1} binary;                # <--: link left

              var wr{0..Nr+1,0..Nc+1} binary;                # -->: link right

              var wu{0..Nr+1,0..Nc+1} binary;               #  A : link upper

              var wd{0..Nr+1,0..Nc+1} binary;                #  V : link lower

s.t. INI4{r in 0..Nr+1,c in 0..Nc+1: r=0||r=Nr+1||c=0||c=Nc+1}: # outer

                             wl[r,c]+wr[r,c]+wu[r,c]+wd[r,c]=0;

s.t. WHT01{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}:                           # center  --> 0 arrow

                             wl[r,c]+wr[r,c]+wu[r,c]+wd[r,c]=0;

s.t. LLL02{r in 1..Nr,c in 1..Nc:ini[r,c]=0}:                             # not center-->1 arrow

                             wr[r,c]+wl[r,c]+wu[r,c]+wd[r,c]=w[r,c];

s.t. LLL03{r in 1..Nr,c in 1..Nc}: wr[r,c]<=w[r,c+1];   # w<--

s.t. LLL04{r in 1..Nr,c in 1..Nc}: wl[r,c]<=w[r,c-1];    #  -->w

s.t. LLL05{r in 1..Nr,c in 1..Nc}: wu[r,c]<=w[r-1,c];   # w  V

s.t. LLL06{r in 1..Nr,c in 1..Nc}: wd[r,c]<=w[r+1,c];   # A  w

 

 ここのコードと、後続のコードで個別の白マスの個数を数え上げる。まず、木構造の構築。

 変数wl[,]wr[,],wu[,]wd[,]は、個別のマスからでていく矢印。マス[r,c]から左にでる矢印がwl[r,c]。以下も同様で、説明は略す。

 制約INI4は、拡大した部分のマスからでていく矢印はない、というもの。この制約、別の制約との絡みで不要とも思えるが、無いと遅くなるので残してある。

 制約WHT01は、数字指定のマスから出て行く矢印はない、というもの。無くてもいかとも思うが、いれてある。

 制約LLL02は、数字指定のマス以外であれば、白マスであれば1つ出て行く矢印があるというもの。黒マスであれば出て行く矢印はない。

 制約LLL03などは、矢印が出て行った先のマスは白色というもの。制約LLL03は右矢印。制約LLL04は左矢印。以下2つは略。

 

## ------------ count white cells in the same islands ------------------

              var nl{0..Nr+1,0..Nc+1} >=0;

              var nr{0..Nr+1,0..Nc+1} >=0;

              var nu{0..Nr+1,0..Nc+1} >=0;

              var nd{0..Nr+1,0..Nc+1} >=0;

s.t. MEMBl{r in 0..Nr+1,c in 0..Nc+1}: nl[r,c]<=N*wl[r,c];

s.t. MEMBr{r in 0..Nr+1,c in 0..Nc+1}: nr[r,c]<=N*wr[r,c];

s.t. MEMBu{r in 0..Nr+1,c in 0..Nc+1}: nu[r,c]<=N*wu[r,c];

s.t. MEMBd{r in 0..Nr+1,c in 0..Nc+1}: nd[r,c]<=N*wd[r,c];

 

s.t. WHT50{r in 1..Nr-0,c in 1..Nc-0}:                    # gather numbers from sons

              wban[r,c]=w[r,c]+nr[r,c-1]+nl[r,c+1]+nd[r-1,c]+nu[r+1,c];

s.t. WHT51{r in 1..Nr-0,c in 1..Nc:ini[r,c]=0}:           # allocate number to one arrow

              wban[r,c]=       nl[r,c]+nr[r,c]+nu[r,c]+nd[r,c];

 

 変数nl[,],nr[,],nu[,],nd[,]は、親に渡すための個数。実数である点に注意。整数でもいいのだが、そうすると極端に遅くなるので、実数にしてある。

 制約MEMBlは、wl[r,c]が0ならnl[r,c]も0であるとするもの。以下の3つの制約も同様。

 制約WHT50は、子からの個数を加算し、さらに自分の1個を加算し、自分の個数とするもの。自分が白でないなら、w[r,c]は0だし、子もいないのでwban[r,c]=0となる。

 制約WHT51は、自分の数えた個数を親に渡すもの。出て行く矢印wl[r,c]などは1つだけ(0でなく)1であり、一人の親に、数えた個数を全て渡すことになる。

 なお、ここのやり方、巡回セールスマン問題でのやり方(Andrew Makhorin)を参考にしている。

 以上で、白マスの部分の説明が終わった。残りは、黒マス。

 

## +++++++++++ the Black Current  ++++++++++++++++++++++++++

              var blr{1..Nr,1..Nc} >=0,<=1;                                  # black link,left right

              var bud{1..Nr,1..Nc} >=0,<=1;                                 # black link, up down

s.t. BLACK1{r in 1..Nr,c in 1..Nc}: blr[r,c]<=b[r,c];

s.t. BLACK2{r in 1..Nr,c in 1..Nc}: blr[r,c]<=1-w[r,c-1];

s.t. BLACK3{r in 1..Nr,c in 1..Nc}: bud[r,c]<=b[r,c];

s.t. BLACK4{r in 1..Nr,c in 1..Nc}: bud[r,c]<=1-w[r-1,c];

 

 変数blr[,]は、自分のマスと左のマスが共に黒かどうかを示す。変数bud[,]は縦方向。

 制約BLACK1,2は、blr[r,c]が1なら、自分は黒で、左のマスは白でない、というもの。

 制約BLACK3,4は、bud[r,c]が1なら、自分は黒で、上のマスは白でない、というもの。

 

 次に、黒マスの数を数える部分が必要だが、その提示は省略。勿論、付録のコードには載せてある。

 

さて、最後に走行とダンプの指定部分。

 

## ++++++++++ solve and dump ++++++++++++++++++++++++++++

              var numB;

s.t. NAMB0: numB=Nr*Nc-sum{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}ini[r,c];

s.t. BHT52{r in 1..Nr,c in 1..Nc:CB&& ini[r,c]=-1}: bban[r,c]=numB;        

              var nnn >=0;

s.t. NNN: nnn=sum{r in 1..Nr,c in 1..Nc}(blr[r,c]+bud[r,c]);

maximize MMM: nnn;

solve;

 

printf("--ini--\n");

for {r in 1..Nr}{ 

  for{c in 1..Nc}{ printf"%3.0f",ini[r,c];} printf("\n");

}

printf "--wban---white=%3.0f--\n",Nr*Nc-numB;

for {r in 1..Nr}{ 

  for{c in 1..Nc}{ printf"%3.0f",wban[r,c];} printf("\n");

}

printf "--bban---black=%3.0f--\n",numB;

for {r in 1..Nr}{ 

  for{c in 1..Nc}{ printf"%3.0f",bban[r,c];} printf("\n");

}

end;

 

変数numBは、黒マスの数を数えるもの。

制約NAMB0は、問題が定める黒マス数を求めている。問題定義配列中の正整数を全て加算すれば、それは白マス数である。それを全マス数Nr*Ncから引けば黒マス数となる。

制約BHT52は、木構造に基づいて数えあげた黒マス数(ini[r,c]=-1の位置のbban[r,c])numBと同じ数になるという制約。丁度、分断禁止制約となる。

次の制約NNNは、黒マス間の隣接個数(blr[,],bud[,])を数え上げて、それを最大化させる。これでほぼ「分断禁止制約」となるはずなのだが、今回は上の分断禁止制約BHT52により、補助的制約に格下げとなった。ただし高速化には貢献していると思う。

 

そして、求解せよ(solve)。終わるとダンプ(printf)。そういうもの。細かい説明は略。

以上で、コードの基本部分は終わった。全コードを説明した気もするのだが、全コードとはいえない。なぜかというと、遅すぎるからである。速くするコードが必要である。

 

[高速化コード]

 ぬりかべパズル、実際にやってみれば分かるのだが、数字を指定した番号の周りから、特に角や辺に配置された数字の周りから、順次に白マス、黒マスが確定していく。簡単な問題では、全てが確定していく。

 

 

 左の2つは、2つの数字指定マスに挟まれたマスは黒マスになるというもの。右の2つは、2つの数字指定マスが頂点を共有して対角に配置されていれば、2つに隣接したマス2つは黒マスになるというもの。

 辺や角であれば、さらに強い制約となる。

 

 

白マスが成長していくパターンも示しておこう。

 

 

こういうパターンをある程度指定しておけば、それなりの速度で走ってくれることになる。ただし、かなり多量に指定しないとならないので、うんざりしてしまう。付録のコードでは、かなり指定してある。うんざりしつつ作業したもの。

 こういう高速化子コード、古い版では角や辺を狙い撃ちしたコードが多かった。新しい版では、一般的に記述するようにした。適用範囲が拡大し、高速化に寄与できたと思う。ただし、コードサイズは膨らんでいる。

 

[ぬりかべとMathProgの相性など]

 最初、分断禁止制約さえ解決すれば良いと思っていた。これは、[解法コード]節の最後のコード内に定義してある変数nnnを最大化すればいいはずだとにらんだ。走行実験の結果、これでいいとの結論になった。なお、新しい版では、分断禁止制約は追加している。

 ただし、高速化の手当てを入れないと遅すぎた。従って[高速化コード]を大量に追加した。という意味で、ぬりかべとMathProgの相性はよくない。ただこれは、どの言語で書いても同様な気がする。

 

[走行例]

 まず7*7のもの。

 

左から、問題、白いシマの数えあげ、黒いマスの数え上げ。

次に8*8

次に10*10

次の10*10

この問題、白いマスが塊になっていたり、数字指定マスがシマの端ではない位置にあったりする。[問題の記述、データ構造]節で示していた状況である。

最後の10*10.

 

 求解時間などは以下。まず、旧版でのもの。

 

問題

p1 7*7

p2 8*8

p3 10*10

p4 10*10

p5 10*10

時間(sec)

0

0.8

9.5

0.3

300以上

2.5GHz

容量(MB)

1.7

2.3

3.4

3

?

2MB*2キャッシュ

 

次に新しい版でのもの。

 

問題

p1 7*7

p2 8*8

p3 10*10

p4 10*10

p5 10*10

時間(sec)

0

0.1

1.3

0.5

1.8

2.2GHz

容量(MB)

3.3

4.1

6

5.8

6.6

i7-2670QM

 

問題P410*10問題なのだが、高速に求解できている。というのは、辺に数字が多く置かれているからである(解答図参照)

問題p5は、古い版では5分経っても終わらなかった。しかし、新しい版では即座に求解できている。新しい版、黒マスの生成や白マスの生成、それなりに素直に指定できたからであろう。

 

[MathProg求解コード]

 

        param CB:=1;                                             #

## +++++++++++ white islands in the Balck Current ++++++++++++++++

        param Nr integer;                                        # number of rows

        param Nc integer;                                        # number of columns

        param N:=Nr*Nc;

        param ini{1..Nr,1..Nc} integer;                          # board to define problem

        param SN:=0.001;

        var w{0..Nr+1,0..Nc+1} binary;                           # 1/0: white/not

        var b{0..Nr+1,0..Nc+1} binary;                           # 1/0: black/not

## ------------ initialization ---------------------------

        var wban{0..Nr+1,0..Nc+1} >=0;                           # work board for counting white cells

s.t. INI0{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}: wban[r,c]=ini[r,c];#

s.t. INI2{r in 0..Nr+1,c in 0..Nc+1:                             # outers cell: black and white

                r=0||r=Nr+1||c=0||c=Nc+1}: b[r,c]+w[r,c]=2;

s.t. BLACK0{r in 1..Nr,c in 1..Nc}: b[r,c]=1-w[r,c];             # black is "not white" on the board.

## ------------ inhibit vortexes caused by black current ---------------

s.t. VORT{r in 2..Nr,c in 2..Nc}: b[r-1,c-1]+b[r-1,c]+b[r,c-1]+b[r,c]<=3;        # bB

## ------------ optimizations for cutting redundant search ----------------

s.t. OPT1w{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}: w[r,c]=1; # center cell: white    # n

s.t. OPT1b{r in 1..Nr,c in 1..Nc:ini[r,c] =1}: b[r-1,c]+b[r,c+1]+b[r+1,c]+b[r,c-1]=4;     # 1

s.t. OPT2c{r in 1..Nr,c in 1..Nc:ini[r,c]>=2}: b[r-1,c]+b[r,c+1]+b[r+1,c]+b[r,c-1]<=3;    # 2

s.t. OPT2d{r in 1..Nr,c in 2..Nc:ini[r,c] =2||ini[r,c-1] =2}:                             # 2W

        6*w[r,c]+6*w[r,c-1]<=6+b[r,c-2]+b[r-1,c-1]+b[r-1,c]+b[r,c+1]+b[r+1,c]+b[r+1,c-1];

s.t. OPT2e{r in 2..Nr,c in 1..Nc:ini[r,c] =2||ini[r-1,c] =2}:                               # 2

        6*w[r,c]+6*w[r-1,c]<=6+b[r-2,c]+b[r-1,c-1]+b[r,c-1]+b[r+1,c]+b[r,c+1]+b[r-1,c+1]; # W

 

s.t. OPT3f{r in 1..Nr,c in 2..Nc:ini[r,c]>=3||ini[r,c-1]>=3}:                             # 3W

        b[r,c-2]+b[r-1,c-1]+b[r-1,c]+b[r,c+1]+b[r+1,c]+b[r+1,c-1]<=5;

s.t. OPT3g{r in 2..Nr,c in 1..Nc:ini[r,c]>=3||ini[r-1,c]>=3}:                             # 3

        b[r-2,c]+b[r-1,c-1]+b[r,c-1]+b[r+1,c]+b[r,c+1]+b[r-1,c+1]<=5;                    # W

s.t. OPT3d{r in 1..Nr,c in 3..Nc:ini[r,c] =3||ini[r,c-1] =3||ini[r,c-2]=3}:

        8*w[r,c]+8*w[r,c-1]+8*w[r,c-2]<=                                         # 3wW

        16+b[r,c-3]+b[r-1,c-2]+b[r-1,c-1]+b[r-1,c]+b[r,c+1]+b[r+1,c]+b[r+1,c-1]+b[r+1,c-2];

s.t. OPT4d{r in 1..Nr,c in 3..Nc:ini[r,c]>=4||ini[r,c-1]>=4||ini[r,c-2]>=4}:              # 4wW

        b[r,c-3]+b[r-1,c-2]+b[r-1,c-1]+b[r-1,c]+b[r,c+1]+b[r+1,c]+b[r+1,c-1]+b[r+1,c-2]<=7;

s.t. OPT00{r in 1..Nr,c in 2..Nc-1:ini[r,c-1]>=1&&ini[r,c+1]>=1}: b[r,c]=1;              # n.m

s.t. OPT11{r in 2..Nr-1,c in 1..Nc:ini[r-1,c]>=1&&ini[r+1,c]>=1}: b[r,c]=1;               # n.m  vertical

s.t. OPT2{r in 2..Nr,c in 2..Nc:                                                 # n

                ini[r-1,c-1]>=1&&ini[r,c]>=1}:                                           #  M

                b[r,c-1]+b[r-1,c]=2;                                                    

s.t. OPT3{r in 2..Nr,c in 2..Nc:                                                 #  n

                ini[r,c-1]>=1&&ini[r-1,c]>=1}:                                           # m.

                b[r-1,c-1]+b[r,c]=2;

 

## ------------ id numbers for white cells ------------------------------------

        var wid{0..Nr+1,0..Nc+1} >=0;            # id number for white cells

s.t. WID00{r in 0..Nr+1,c in 0..Nc+1: r=0||r=Nr+1||c=0||c=Nc+1}: wid[r,c]=0;     # outer cells have not id number

s.t. WID01{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}:   wid[r,c]=Nc*(r-1)+c;             # center cells have id number

s.t. WID02{r in 1..Nr,c in 1..Nc}: wid[r,c]<=N*w[r,c];                           # black cells have not id number

 

s.t. WID038{r in 1..Nr,c in 2..Nc}: wid[r,c-1]<=N*(b[r,c]+b[r,c-1])+wid[r,c];             # wW-->wid=wid

s.t. WID039{r in 1..Nr,c in 2..Nc}: wid[r,c  ]<=N*(b[r,c]+b[r,c-1])+wid[r,c-1]; 

s.t. WID040{r in 2..Nr,c in 1..Nc}: wid[r-1,c]<=N*(b[r,c]+b[r-1,c])+wid[r,c];             # w --> wid=wid

s.t. WID041{r in 2..Nr,c in 1..Nc}: wid[r,  c]<=N*(b[r,c]+b[r-1,c])+wid[r-1,c];           # W

 

s.t. OPTd0{r in 2..Nr,c in 2..Nc}: w[r-1,c-1]+w[r,c]+SN*(wid[r-1,c-1]-wid[r,c])<=2+b[r-1,c];      # wid

s.t. OPTd1{r in 2..Nr,c in 2..Nc}: w[r-1,c-1]+w[r,c]+SN*(wid[r,c]-wid[r-1,c-1])<=2+b[r-1,c];      #    Wid

 

s.t. OPTd2{r in 2..Nr,c in 2..Nc}: w[r-1,c-1]+w[r,c]+SN*(wid[r-1,c-1]-wid[r,c])<=2+b[r,c-1];      #    wid

s.t. OPTd3{r in 2..Nr,c in 2..Nc}: w[r-1,c-1]+w[r,c]+SN*(wid[r,c]-wid[r-1,c-1])<=2+b[r,c-1];      # wid .

 

s.t. OPTd4{r in 1..Nr,c in 2..Nc-1}: w[r,c-1]+w[r,c+1]+SN*(wid[r,c-1]-wid[r,c+1])<=2+b[r,c]; # wid . wid

s.t. OPTd5{r in 2..Nr-1,c in 1..Nc}: w[r-1,c]+w[r+1,c]+SN*(wid[r-1,c]-wid[r+1,c])<=2+b[r,c]; # same,vertical

 

## ------------ arrows for counting white cells ---------------------------

        var wl{0..Nr+1,0..Nc+1} binary;                  # <--: link left

        var wr{0..Nr+1,0..Nc+1} binary;                  # -->: link right

        var wu{0..Nr+1,0..Nc+1} binary;                  #  ^ : link upper

        var wd{0..Nr+1,0..Nc+1} binary;                  #  v : link lower

s.t. INI4{r in 0..Nr+1,c in 0..Nc+1: r=0||r=Nr+1||c=0||c=Nc+1}: # outer cells

                wl[r,c]+wr[r,c]+wu[r,c]+wd[r,c]=0;

s.t. WHT01{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}:           # center  --> 0 arrow

                wl[r,c]+wr[r,c]+wu[r,c]+wd[r,c]=0;

s.t. LLL02{r in 1..Nr,c in 1..Nc:ini[r,c]=0}:            # not center-->1 arrow

                wl[r,c]+wr[r,c]+wu[r,c]+wd[r,c]=w[r,c];

s.t. LLL03{r in 0..Nr+1,c in 0..Nc+0}: wr[r,c]<=1-b[r,c+1];      # -->w

s.t. LLL04{r in 0..Nr+1,c in 1..Nc+0}: wl[r,c]<=1-b[r,c-1];      # w<--

s.t. LLL05{r in 1..Nr+1,c in 0..Nc+1}: wu[r,c]<=1-b[r-1,c];      # w  v

s.t. LLL06{r in 0..Nr+0,c in 0..Nc+1}: wd[r,c]<=1-b[r+1,c];      # ^  w

## ------------ count white cells in the same islands ------------------

        param NN:=max{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}ini[r,c];

        var nl{0..Nr+1,0..Nc+1} >=0,<=NN;

        var nr{0..Nr+1,0..Nc+1} >=0,<=NN;

        var nu{0..Nr+1,0..Nc+1} >=0,<=NN;

        var nd{0..Nr+1,0..Nc+1} >=0,<=NN;

s.t. MEMBl{r in 0..Nr+1,c in 0..Nc+1}: nl[r,c]<=NN*wl[r,c];

s.t. MEMBr{r in 0..Nr+1,c in 0..Nc+1}: nr[r,c]<=NN*wr[r,c];

s.t. MEMBu{r in 0..Nr+1,c in 0..Nc+1}: nu[r,c]<=NN*wu[r,c];

s.t. MEMBd{r in 0..Nr+1,c in 0..Nc+1}: nd[r,c]<=NN*wd[r,c];

 

s.t. WHT50{r in 1..Nr-0,c in 1..Nc-0}:           # gather numbers from sons

        wban[r,c]=w[r,c]+nr[r,c-1]+nl[r,c+1]+nd[r-1,c]+nu[r+1,c];

s.t. WHT51{r in 1..Nr-0,c in 1..Nc:ini[r,c]=0}:  # allocate number to one arrow

        wban[r,c]=       nl[r,c]+nr[r,c]+nu[r,c]+nd[r,c];

 

## +++++++++++ the Black Current  ++++++++++++++++++++++++++

        var blr{1..Nr,1..Nc} >=0,<=1;                    # black link,left right

        var bud{1..Nr,1..Nc} >=0,<=1;                    # black link, up down

s.t. BLACK1{r in 1..Nr,c in 1..Nc}: blr[r,c]<=b[r,c];

s.t. BLACK2{r in 1..Nr,c in 1..Nc}: blr[r,c]<=1-w[r,c-1];

s.t. BLACK3{r in 1..Nr,c in 1..Nc}: bud[r,c]<=b[r,c];

s.t. BLACK4{r in 1..Nr,c in 1..Nc}: bud[r,c]<=1-w[r-1,c];

 

## ------------ arrows for counting black cells ---------------------------

        var bl{0..Nr+1,0..Nc+1} binary;                  # <--: link left

        var br{0..Nr+1,0..Nc+1} binary;                  # -->: link right

        var bu{0..Nr+1,0..Nc+1} binary;                  #  ^ : link upper

        var bd{0..Nr+1,0..Nc+1} binary;                  #  v : link lower

s.t. BNI4{r in 0..Nr+1,c in 0..Nc+1:CB&& (r=0||r=Nr+1||c=0||c=Nc+1)}: # outer cells

                bl[r,c]+br[r,c]+bu[r,c]+bd[r,c]=0;

s.t. BHT01{r in 1..Nr,c in 1..Nc:CB&& ini[r,c]=-1}:              # center  --> 0 arrow

                bl[r,c]+br[r,c]+bu[r,c]+bd[r,c]=0;

s.t. BLL02{r in 1..Nr,c in 1..Nc:CB&& ini[r,c]=0}:               # not center-->1 arrow

                bl[r,c]+br[r,c]+bu[r,c]+bd[r,c]=b[r,c];

s.t. BLL03{r in 0..Nr+1,c in 0..Nc+0:CB}: br[r,c]<=1-w[r,c+1];   # -->w

s.t. BLL04{r in 0..Nr+1,c in 1..Nc+0:CB}: bl[r,c]<=1-w[r,c-1];   # w<--

s.t. BLL05{r in 1..Nr+1,c in 0..Nc+1:CB}: bu[r,c]<=1-w[r-1,c];   # w  v

s.t. BLL06{r in 0..Nr+0,c in 0..Nc+1:CB}: bd[r,c]<=1-w[r+1,c];   # ^  w

## ------------ count black cells in the Black Current ------------------

        var bban{0..Nr+1,0..Nc+1} >=0;           # work board for counting black cells

        var nbl{0..Nr+1,0..Nc+1} >=0,<=Nr*Nc;    # transported number with bl

        var nbr{0..Nr+1,0..Nc+1} >=0,<=Nr*Nc;    # transportrd number with br

        var nbu{0..Nr+1,0..Nc+1} >=0,<=Nr*Nc;    # transportrd number with bu

        var nbd{0..Nr+1,0..Nc+1} >=0,<=Nr*Nc;    # transportrd number with bd

s.t. BEMBl{r in 0..Nr+1,c in 0..Nc+1:CB}: nbl[r,c]<=Nr*Nc*bl[r,c];

s.t. BEMBr{r in 0..Nr+1,c in 0..Nc+1:CB}: nbr[r,c]<=Nr*Nc*br[r,c];

s.t. BEMBu{r in 0..Nr+1,c in 0..Nc+1:CB}: nbu[r,c]<=Nr*Nc*bu[r,c];

s.t. BEMBd{r in 0..Nr+1,c in 0..Nc+1:CB}: nbd[r,c]<=Nr*Nc*bd[r,c];

s.t. BHT50{r in 1..Nr-0,c in 1..Nc-0:CB}:                # gather numbers from sons

        bban[r,c]=b[r,c]+nbr[r,c-1]+nbl[r,c+1]+nbd[r-1,c]+nbu[r+1,c];

s.t. BHT51{r in 1..Nr-0,c in 1..Nc:CB&& ini[r,c]=0}:     # allocate number to one arrow

        bban[r,c]=       nbl[r,c]+nbr[r,c]+nbu[r,c]+nbd[r,c];

 

## +++++++++++ inhibit black pool ++++++++++++++++++++++++++

s.t. GenB10{r in 1..Nr,c in 1..Nc}:

        w[r,c-1]+w[r-1,c]+w[r,c+1]+w[r+1,c]+b[r,c]<=4;                   # b

s.t. GenB20{r in 1..Nr,c in 1..Nc-1}:                                    # 

        w[r,c-1]+w[r-1,c]+w[r-1,c+1]+w[r,c+2]+w[r+1,c+1]+w[r+1,c]+       # Bb

        0.5*(b[r,c]+b[r,c+1])<=6;                                        # 

s.t. GenB21{r in 1..Nr-1,c in 1..Nc}:                                    #

        w[r-1,c]+w[r,c+1]+w[r+1,c+1]+w[r+2,c]+w[r+1,c-1]+w[r,c-1]+       # B

        0.5*(b[r,c]+b[r+1,c])<=6;                                        # b

s.t. GenB30{r in 1..Nr,c in 1..Nc-2}:                           

        w[r-1,c]+w[r-1,c+1]+w[r-1,c+2]+w[r,c+3]+w[r+1,c+2]+w[r+1,c+1]+   # Bbb

        w[r+1,c]+w[r,c-1]+0.3*(b[r,c]+b[r,c+1]+b[r,c+2])<=8;            

s.t. GenB31{r in 1..Nr-2,c in 1..Nc}:                                    # B    

        w[r-1,c]+w[r,c+1]+w[r+1,c+1]+w[r+2,c+1]+w[r+3,c]+w[r+2,c-1]+     # b

        w[r+1,c-1]+w[r,c-1]+0.3*(b[r,c]+b[r+1,c]+b[r+1,c])<=8;           # b

s.t. GenB32{r in 2..Nr,c in 2..Nc}: w[r-1,c-1]+w[r-2,c]+ #  

        w[r-1,c+1]+w[r,c+1]+w[r+1,c]+w[r+1,c-1]+w[r,c-2]+        #  b

        0.3*(b[r-1,c]+b[r,c-1]+b[r,c])<=7;                       # bB

s.t. GenB33{r in 2..Nr,c in 2..Nc}: w[r-2,c-1]+w[r-1,c]+ #  

        w[r,c+1]+w[r+1,c]+w[r+1,c-1]+w[r,c-2]+w[r-1,c-2]+        # b

        0.3*(b[r-1,c-1]+b[r,c-1]+b[r,c])<=7;                     # bB

s.t. GenB34{r in 2..Nr,c in 2..Nc}: w[r-2,c-1]+w[r-2,c]+ #  

        w[r-1,c+1]+w[r,c]+w[r+1,c-1]+w[r,c-2]+w[r-1,c-2]+        # bb

        0.3*(b[r-1,c-1]+b[r,c-1]+b[r-1,c])<=7;                   # b.

s.t. GenB35{r in 2..Nr,c in 2..Nc}: w[r-2,c-1]+w[r-2,c]+ #  

        w[r-1,c+1]+w[r,c+1]+w[r+1,c]+w[r,c-1]+w[r-1,c-2]+        # bb

        0.3*(b[r-1,c-1]+b[r,c]+b[r-1,c])<=7;                     #  B

s.t. GenB40{r in 3..Nr-1,c in 3..Nc-1}: w[r-2,c-1]+w[r-2,c]+     #

        w[r-1,c+1]+w[r,c+1]+w[r+1,c]+w[r+1,c-1]+w[r,c-2]+        # bb

        w[r-1,c-2]+0.25*(b[r-1,c-1]+b[r-1,c]+b[r,c-1]+b[r,c])<=8;# bB

 

## ++++++++++ solve and dump ++++++++++++++++++++++++++++

        var numB;

s.t. NAMB0: numB=Nr*Nc-sum{r in 1..Nr,c in 1..Nc:ini[r,c]>=1}ini[r,c];

s.t. BHT52{r in 1..Nr,c in 1..Nc:CB&& ini[r,c]=-1}: bban[r,c]=numB;     

        var nnn >=0;

s.t. NNN: nnn=sum{r in 1..Nr,c in 1..Nc}(blr[r,c]+bud[r,c]);

maximize MMM: nnn;

solve;

 

printf("--ini--\n");

for {r in 1..Nr}{ 

  for{c in 1..Nc}{ printf"%3.0f\t",ini[r,c];} printf("\n");

}

printf "--wban---white=%3.0f--\n",Nr*Nc-numB;

for {r in 1..Nr}{ 

  for{c in 1..Nc}{ printf"%3.0f\t",wban[r,c];} printf("\n");

}

printf "--bban---black=%3.0f--\n",numB;

for {r in 1..Nr}{ 

  for{c in 1..Nc}{ printf"%3.0f\t",bban[r,c];} printf("\n");

}

end;

 

 

[問題定義コード]

 

#p01.txt

data;

         param Nr:=7;

         param Nc:=7;

         param ini:

          1 2 3 4 5 6 7:=

  1      0 3 0 0 0 1 -1

  2      0 0 0 0 0 0 0

  3      2 0 0 1 0 0 0

  4      0 0 0 0 0 0 0

  5      0 1 0 0 2 0 0

  6      0 0 2 0 0 0 0

  7      1 0 0 0 1 0 6

; end;

 

#p02.txt

data;

        param Nr:=8;

        param Nc:=8;

        param ini:

         1 2 3 4 5 6 7 8:=

  1     1 0 0 0 0 2 0 0

  2     0 0 0 1 0 0 0 0

  3     0 0 0 0 0 0 0 3

  4     0 3 0 0 0 0 0 0

  5     0 0 0 0 0 0 0 0

  6     0 0 0 0 0 0 0 0

  7     0 0 0 0 0 4 -1 0

  8     4 0 6 0 0 0 5 0

; end;

 

 

#p03.txt

data;

        param Nr:=10;

        param Nc:=10;

        param ini:

         1 2 3 4 5 6 7 8 9 10:=

  1     0 0 0 0 0 0 0 5 -1 2

  2     3 0 0 0 0 0 0 0 0 0

  3     0 4 0 0 2 0 0 0 0 0

  4     0 0 0 0 0 0 3 0 0 0

  5     0 4 0 0 0 4 0 0 0 0

  6     0 0 0 0 0 0 0 0 0 3

  7     0 0 0 0 0 0 0 0 0 0

  8     0 0 0 0 0 0 0 0 0 0

  9     0 3 0 0 3 0 0 0 0 0

  10    0 0 1 0 0 1 0 3 0 3

; end;

 

#p04.txt

data;

        param Nr:=10;

        param Nc:=10;

        param ini:

         1 2 3 4 5 6 7 8 9 10:=

  1     6 0 2 0 3 0 0 0 0 3

  2     0 0 0 0 0 0 0 0 0 -1

  3     0 0 0 0 0 0 0 0 0 4

  4     0 0 0 0 0 0 0 0 0 0

  5     0 0 0 0 2 0 0 0 0 2

  6     3 0 0 0 0 5 0 0 0 0

  7     0 0 0 0 0 0 0 0 0 0

  8     3 0 0 0 0 0 0 0 0 0

  9     0 0 0 0 0 0 0 0 0 0

  10    4 0 0 0 0 5 0 4 0 1

; end;

 

#p05.txt

data;

        param Nr:=10;

        param Nc:=10;

        param ini:

         1 2 3 4 5 6 7 8 9 10:=

  1     0 3 0 0 0 0 4 0 0 0

  2     0 0 0 0 0 6 0 0 0 0

  3     0 0 0 0 0 0 0 2 0 0

  4     0 0 0 0 0 0 3 0 0 0

  5     0 0 0 0 0 0 0 0 2 0

  6     0 4 0 0 0 0 0 3 0 -1

  7     0 0 0 0 0 0 0 0 0 1

  8     0 10 0 0 0 0 0 0 3 0

  9     0 0 0 0 0 0 0 0 0 0

  10    0 0 3 0 0 0 0 0 0 2

; end;

 

 

トップページへ

 

Ads by TOK2