GLPKでの「ひとりにしてくれ」の攻略法(解法)

キューブ王

パズルの問題を解くプログラムに挑戦している。今回、「ひとりにしてくれ」に挑戦。GLPKMathProgで書いた。なお前回、スリリンに挑戦した。旨く攻略できている。

この「ひとりにしてくれ」、領域が分断されてはならないという制約条件がある。分断系の制約条件、整数計画法で記述するのは面倒だと言われている。しかし、スリリンのループ拒否コードと同様なやり方が可能なはずで挑戦してみたもの。簡単に攻略できた。その紹介。

 

[ひとりにしてくれのルール]

 このパズルのルールをwikipediaに従い説明しておく。以下の左図が「ひとりにしてくれ」の盤面例。5*5の盤面である。

 

 

(1) 盤面に黒マスを配置して、各行・各列に同じ数字が残らないようにする。

(2) 黒マスは縦または横に連続しない。

(3) また、黒マスで盤面が分断されてはいけない。

「ひとりにしてくれ」はこの2,3番目のルール(連黒分断禁と呼ばれる基本的な黒マスルール)をはじめて取り入れたパズルでもある。

 

[直ぐ分かる攻略法]

人間が手作業で解くとすれば、以下の攻略法を直ぐに思いつく。

(a) 同じ数が3つ並んでいると、両側は黒マスになる。

(b) 同じ数が2つ並んでいると、片方が白マスになる。 さらに、同じ列にその数字はもう入らない。

これらを入口にして、黒マスになるマス・白マスになるマスを決めていくことで解答が得られる。例題では上部の1の集まりがこれに当たる。

中級以上の問題になると、「どちらを黒マスにしてもこのマスは白マスになる」という議論を経る事や、「全体を分断する黒マスの並びを回避するための白マス」を探す事が多くなる。

 

 以上、wikipediaの説明をほぼ転載。右の図は解答例である。解は1つという制約があるはずで、右の図は解等例ではなく唯一解という方が正しいか。

 

 

 

[整数計画法での解法(制約の記述法)]

 さて、この問題の制約条件を記述するには、その前に、マスやノードを特定するための名称系や座標系を定める必要がある。例えば、以下でいい。一般にはN*Nの盤面なのだが、2*2で示す。2*2の上下左右に余分にマスを追加して4*4になっている。

 

 

本来の盤面は、黄色いマス4つ(2*2)である。ノード(小さい四角)は植木算で9(3*3)である。個別の片に付けた1,1等の数値の組みは、個別片を識別する座標である。ノードは、そこを通過する斜線を識別するためのもの。dur, ddl, ddr, dulの4種ある。可能な方向4種に対応する。

 次に本論の制約条件。これは1つを除いて非常に簡単。1つとは、[ひとりにしてくれのルール](3)に対応するもの。つまり、黒マスで盤面が分断されてはいけないというもの。これは後で説明する。

 

## ----------- alone -------------------

           param N;

           param M:= 2*N;

           var bbb{0..N+1,0..N+1,0..N} ,binary;   # work board

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

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

s.t. INI0{ r in 0..N+1 by N+1, c in 0..N+1}:                 bbb[r,c,0] =1;      # outer

s.t. INI1{ r in 0..N+1 by N+1, c in 0..N+1, n in 1..N}: bbb[r,c,n] =0;          # outer

s.t. INI2{ r in 0..N+1, c in 0..N+1 by N+1}:                 bbb[r,c,0] =1;      # outer

s.t. INI3{ r in 0..N+1, c in 0..N+1 by N+1, n in 1..N}: bbb[r,c,n] =0;          # outer

s.t. INI5{ r in 1..N, c in 1..N, n in 1..N: n= ini[r,c]}: bbb[r,c,n] <=1; # copy

s.t. INI6{ r in 1..N, c in 1..N, n in 1..N: n!=ini[r,c]}: bbb[r,c,n]  =0; # copy

 

param Nは問題の規模。N*N行列での問題ということを示す。問題を解く際に使う3次元配列がbbb[,,]

param ini[,]は、問題定義配列。上の例でいえば、黄色い部分に相当。個別のマスの数字を指定(問題定義)

制約INI0-INI3は、本来の盤面の外側に0を設定するもの。ここでの0とは、黒マスにするという意味である。

制約INI5,INI6は、本来の盤面内のマス[r,c]ini[r,c]または0(黒マス)を設定するもの。

 

## ------------ nearly equal to latin square -------------

s.t. LAT0{ r in 1..N  , c in 1..N  }: sum{ n in 0..N} bbb[r,c,n]  =1;

s.t. LAT1{ r in 1..N  , n in 1..N  }: sum{ c in 1..N} bbb[r,c,n] <=1;

s.t. LAT2{ c in 1..N  , n in 1..N  }: sum{ r in 1..N} bbb[r,c,n] <=1;

s.t. LAT3{ r in 1..N  , c in 1..N-1}: bbb[r,c,0]+bbb[r,c+1,0] <=1;

s.t. LAT4{ r in 1..N-1, c in 1..N  }: bbb[r,c,0]+bbb[r+1,c,0] <=1;

 

制約LAT0は、個別のマス[r,c]に1つの数字nを設定するもの。制約INI5,INI6を考えれば、nini[r,c]か0である。

制約LAT1,LAT2は、[ひとりにしてくれのルール](1)に対応。各行・各列に同じ数字が残らないようにするもの。

制約LAT3,LAT4は、[ひとりにしてくれのルール](2)に対応。黒マスは縦横に連続しないようにするもの。

 

[分断禁止制約の指定法]

 分断禁止の制約、どう指定するか、最初分からなかった。しかし、「分断する斜線系列の禁止」と考えなおせば、ほぼ、「輪となった斜線系列の禁止」と同じとなる。であるなら、巡回セールスマン問題の解法と同様なもので解決できることになる。

 

 

 左端の図では、黒マスは盤面を分断していない。右の2つは分断している。中の図、周辺の黒()マスを含めて考えると輪ができている。結局、輪を拒否すればいいことになる。

 

 

 周辺の灰色マスに全て番号1を与える。また、有向斜線で黒マスを結線していく。辿っていく順に番号を振る。全ての有向斜線で到着先の番号が出発元の番号よりも大きければOK、そうでなければNG。そういう制約でいい。この制約、輪がなければ満たすことができ、輪があれば満たすことができないからである。

 なお、個別のマスについて、入ってくる有効斜線は1つ以下という制約も付いている。この意味を示しておく。

 

 

 左図の矢印群までは制約条件を満たしている。最後の1つの斜め矢印、左上を向いても、右下を向いても、行った先のマスへ2つの斜め矢印が入ることになり制約条件を満たさない。つまりこの図の場合、マスに番号を与えるまでもなく、黒マス群がループを形成していて、従って盤面を分断していることが分かることになる。

 一方、入ってくる有効斜線が1つでなくてもいいとすれば、右図が可能となる。結線群がループを形成し、しかも番号を振ることもできている。このループは本来、拒否されるべきものであった。

 

[分断禁止制約のコード]

コードを示しておこう。

 

## ------------ no separation(directed graph) ----------------------

           var ddl{1..N+1,1..N+1} binary;          # down,left

           var ddr{1..N+1,1..N+1} binary;          # down,right

           var dul{1..N+1,1..N+1} binary;          # up,left

           var dur{1..N+1,1..N+1} binary;          # up,right

s.t. XX0{ r in 1..N+1, c in 1..N+1}: dur[r,c]+ddl[r,c]<=1; # 1 direction

s.t. XX1{ r in 1..N+1, c in 1..N+1}: ddr[r,c]+dul[r,c]<=1; # 1 direction

s.t. XX2{ r in 1..N+0, c in 1..N+0}:

                       ddr[r,c]+dur[r+1,c]+dul[r+1,c+1]+ddl[r,c+1]<=1;          # 1 arrival

s.t. SEP0{ c in 2..N-1}: bbb[1,c,0]=ddr[1,c];                          # outer

s.t. SEP1{ c in 2..N-1}: bbb[N,c,0]=dur[N+1,c];                                  # outer

s.t. SEP2{ r in 1..N-0}: bbb[r,1,0]=dur[r+1,1];                         # outer

s.t. SEP3{ r in 1..N-0}: bbb[r,N,0]=dul[r+1,N+1];                     # outer

s.t. SEP8{ r in 1..N-1, c in 1..N}:                            # arc

                       bbb[r,c,0]+bbb[r+1,c-1,0]<=1+dur[r+1,c]+ddl[r+1,c];

s.t. SEP9{ r in 1..N-1, c in 1..N}:                            # arc

                       bbb[r,c,0]+bbb[r+1,c+1,0]<=1+dul[r+1,c+1]+ddr[r+1,c+1];

 

変数が4つある。ddl, ddr, dul, durである。有向斜線の可能な方向に対応する。

制約XX0は、右上に向かう斜線と左下に向かう斜線の2つがあってはならないというもの。

制約XX1は、右下に向かう斜線と左上に向かう斜線の2つがあってはならないというもの。

制約XX2は、個別のマスに入る斜線は1つだけというもの。

制約SEP0,1,2,3は、本来の盤面(黄色部分)に外から入る斜線を設定するもの(上の図)SEP0は上から、SEP1は下から、SEP2は左から、SEP3は右からである。

制約SEP8,9は、盤面の内側の制約を示す。1つの頂点を共有した2つのマスが共に黒マスなら、2つのマスを結ぶ斜線があるというもの。

 

最後に、盤面内で(辺に接しないで)ループを構成する部分を撥ねつける必要があるのだが、この説明は省略。コードは、このページの最後尾に載っている。

 

 [求解時間など]

 実は、ニコリの「ひとりにしてくれ」本屋に買いにいったのだが、置いてなかったので、買えなかった。それで、インターネットでお試し問題を調べてみた。8*8問題が5つ、次に12*12問題がいくつか、次に17*17問題がいくつか用意されていた。最初の7個をやってみた。全て、正解した(正解するまでデバッグした)

 

問題

p01

p02

p03

p04

p05

p06

p07

sec

0.1

0.0

0.2

0.0

0.1

0.4

2.8

サイズ

2.3M

2.3M

2.4M

2.3M

2.4M

5.2M

5.2M

 

このパズル、予想通り、コンピュータとの相性がいいようである。というか、場合の数が少ないだけか。

 なお、盤面内で(辺に接しないで)ループを構成する部分を撥ねつけなくても、かなりの割合で正解できるように思った。特に8*8問題の場合。

 ループ拒否のコード、巡回セールスマン問題の2番目に相当するコード(スリリン参照)に置き換えれば、12*12問題(P06,P07)も格段に速くなると思う。

 

 そういうことで、コードの微調整を試みた。

 まず、巡回セールスマン問題の2番目相当のコードにトライした、旨くいかなかった。巡回セールスマン問題では、結線が三つ又になることはない。しかし、「ひとりで」では三つ又になりうる。その違いからだろう、旨くいかなかった。

 つぎに、ループを構成する部分の撥ねつけコードを不活性にしてみた。試した問題全てで正解した。走行時間が大幅に削減されるはずだったのだが、実際には遅くなった。整数計画法の素人である当方には、何が起こったのかさっぱり分からない。

 最後に、目標関数を吟味した。盤面の中央部分での結線数が少ないほうが、正解に近いはずとにらんだ。ノード4個の結線がループすれば結線数は4個。ループしてなければ3個だろう。ループ拒否と結線数少は近い関係にある。そういう程度の判断であるが、結線数少の目標関数で、規模が大きい問題で2倍近く速くなった気がする。末尾のコードは、この調整を加えたものを載せている。

 

問題

p01

p02

p03

p04

p05

p06

p07

p08

sec

0.1

0.1

0.0

0.0

0.0

0.2

0.9

0.4

サイズ

2.3M

2.3M

2.4M

2.3M

2.4M

5.3M

5.3M

5.4M

 

p08は、12*12問題。以降、17*17問題が3つあるのだが、問題の入力にうんざりして、それは試してない。

 

[後日の分析] 111206

(1)   矢印の方向を逆にすれば、巡回セールスマンの2番目風のループ拒否コードを記述できた。しかし遅かった。

(2)   ループ拒否コード不活性で遅くなる理由は、小さいループがあれば、目標関数値を最小化できず、従ってループ拒否が、目標値到達を早めるからであろう。ただし、単なる推定である。

 

[ひとりにしてくれのコード]

 

以下のような起動法となる。glpsol445は、標準には/usr/local/bin/glpsolだったか。

 

./glpsol445 --dfs -m $1_m.txt -d $2.txt -o out1.txt

 

問題定義のini[,]は別ファイルとした。p01.txtなど。

 

## ----------- alone -------------------

           param N;

           param M:= 2*N;

           var bbb{0..N+1,0..N+1,0..N} ,binary;   # work board

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

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

s.t. INI0{ r in 0..N+1 by N+1, c in 0..N+1}:                 bbb[r,c,0] =1;      # outer

s.t. INI1{ r in 0..N+1 by N+1, c in 0..N+1, n in 1..N}: bbb[r,c,n] =0;          # outer

s.t. INI2{ r in 0..N+1, c in 0..N+1 by N+1}:                 bbb[r,c,0] =1;      # outer

s.t. INI3{ r in 0..N+1, c in 0..N+1 by N+1, n in 1..N}: bbb[r,c,n] =0;          # outer

s.t. INI5{ r in 1..N, c in 1..N, n in 1..N: n= ini[r,c]}: bbb[r,c,n] <=1; # copy

s.t. INI6{ r in 1..N, c in 1..N, n in 1..N: n!=ini[r,c]}: bbb[r,c,n]  =0; # copy

## ------------ nearly equal to latin square -------------

s.t. LAT0{ r in 1..N  , c in 1..N  }: sum{ n in 0..N} bbb[r,c,n]  =1; # 1 value

s.t. LAT1{ r in 1..N  , n in 1..N  }: sum{ c in 1..N} bbb[r,c,n] <=1; # latin

s.t. LAT2{ c in 1..N  , n in 1..N  }: sum{ r in 1..N} bbb[r,c,n] <=1; # latin

s.t. LAT3{ r in 1..N  , c in 1..N-1}: bbb[r,c,0]+bbb[r,c+1,0]<=1; # black black

s.t. LAT4{ r in 1..N-1, c in 1..N  }: bbb[r,c,0]+bbb[r+1,c,0]<=1; # black black

## ------------ no separation(directed graph) ----------------------

           var ddl{1..N+1,1..N+1} binary;          # down,left

           var ddr{1..N+1,1..N+1} binary;          # down,right

           var dul{1..N+1,1..N+1} binary;          # up,left

           var dur{1..N+1,1..N+1} binary;          # up,right

s.t. XX0{ r in 1..N+1, c in 1..N+1}: dur[r,c]+ddl[r,c]<=1; # 1 direction

s.t. XX1{ r in 1..N+1, c in 1..N+1}: ddr[r,c]+dul[r,c]<=1; # 1 direction

s.t. XX2{ r in 1..N+0, c in 1..N+0}:

                       ddr[r,c]+dur[r+1,c]+dul[r+1,c+1]+ddl[r,c+1]<=1;          # 1 arrival

s.t. SEP0{ c in 2..N-1}: bbb[1,c,0]=ddr[1,c];                          # outer

s.t. SEP1{ c in 2..N-1}: bbb[N,c,0]=dur[N+1,c];                                  # outer

s.t. SEP2{ r in 1..N-0}: bbb[r,1,0]=dur[r+1,1];                         # outer

s.t. SEP3{ r in 1..N-0}: bbb[r,N,0]=dul[r+1,N+1];                     # outer

s.t. SEP8{ r in 1..N-1, c in 1..N}:                            # arc

                       bbb[r,c,0]+bbb[r+1,c-1,0]<=1+dur[r+1,c]+ddl[r+1,c];

s.t. SEP9{ r in 1..N-1, c in 1..N}:                            # arc

                       bbb[r,c,0]+bbb[r+1,c+1,0]<=1+dul[r+1,c+1]+ddr[r+1,c+1];

## ------------ no separation(arriving order) ----------------------

           var o1 {0..N+1,0..N+1} >=0, <=M;     # serial number

s.t. NUM0{ r in 0..N+1,                  c in 0..N+1 by N+1}: o1[r,c]=1;       # outer

s.t. NUM1{ r in 0..N+1 by N+1,        c in 0..N+1           }: o1[r,c]=1;     # outer

s.t. NUM4{ r in 1..N, c in 1..N}: o1[r,c]+M*(1-ddr[r,c])    >=o1[r-1,c-1]+1;

s.t. NUM5{ r in 1..N, c in 1..N}: o1[r,c]+M*(1-dur[r+1,c])  >=o1[r+1,c-1]+1;

s.t. NUM6{ r in 1..N, c in 1..N}: o1[r,c]+M*(1-dul[r+1,c+1])>=o1[r+1,c+1]+1;

s.t. NUM7{ r in 1..N, c in 1..N}: o1[r,c]+M*(1-ddl[r,c+1])  >=o1[r-1,c+1]+1;

## ---------- solve and dump ------------------------

           var ban{0..N+1,0..N+1} ,integer,>=0,<=N;

           var obj_value;

s.t. CONV1{ r in 0..N+1, c in 0..N+1}: sum{ n in 0..N}n*bbb[r,c,n] =ban[r,c];

#s.t. MINV: obj_value=sum{ r in 1..N+1, c in 1..N+1}

#s.t. MINV: obj_value=sum{ r in 1..N+0, c in 1..N+0}

s.t. MINV: obj_value=sum{ r in 3..N-1, c in 3..N-1}

                                  (ddr[r,c]+dur[r,c]+dul[r,c]+ddl[r,c]);

 

minimize VALUE: obj_value;

solve;

 

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

  for{c in 0..N+1}{ printf" *" ;} printf"\n";

for {r in 1..N}{ 

  printf" *";

  for{c in 1..N} printf"%2s",if ini[r,c]then ini[r,c] else "@"; printf(" *\n");}

  for{c in 0..N+1}{ printf" *" ;} printf"\n";

 

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

for {r in 0..N+1}{ 

  for{c in 0..N+1}{ printf"%2s",o1[r,c];} printf("\n");}

 

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

  for{c in 0..N+1}{ printf" *" ;} printf"\n";

for {r in 1..N}{ 

  printf" *";

  for{c in 1..N} printf"%2s",if ban[r,c]then ban[r,c] else "@"; printf(" *\n");}

  for{c in 0..N+1}{ printf" *" ;} printf"\n";

end;

 

 

[問題定義ファイル]

 

お試し問題の最初の2つを載せておく。

 

data;      ## ---------- input --------------

param N:= 8;

param ini:           1 2 3 4 5 6 7 8:=

           1          4 5 6 8 3 7 8 1

           2          1 1 1 4 8 6 7 2

           3          5 4 7 6 3 1 2 8

           4         7 4 5 3 2 8 1 6

           5          2 2 2 7 1 4 4 4

           6          8 7 4 5 5 2 6 3

           7          6 7 8 5 5 4 3 2

           8          3 3 3 2 4 6 1 7 ;

end;

 

 

data;      ## ---------- input --------------

param N:= 8;

param ini:           1 2 3 4 5 6 7 8:=

           1          1 1 5 3 3 4 7 8

           2          1 1 6 3 3 8 8 7

           3          7 2 2 4 6 8 8 3

           4         4 2 2 6 1 7 1 5

           5          5 8 7 7 4 3 2 6

           6          8 5 7 7 2 6 4 4

           7          6 6 3 2 5 5 4 4

           8          6 6 4 8 5 5 3 2 ;

end;

 

 

トップページへ

 

Ads by TOK2