GLPKでのペンシルパズル(ナンリンなど)攻略法

キューブ王

0. 導入

パズルの問題を解くプログラムについて考えてみた。普通にはC言語などで書くものだろう。しかしC言語での記述だと、例えばナンバーリンクの解法を書くのはうんざりする。そこでどうしたものか悩んでいた。そういうときに、「数理計画法によるパズル解法」という記事をインターネットで見つけた[]。その記事を読んでみると、なかなか面白かった。その記事では数独やノノグラムの解法作成法が記述してあった。ナンバーリンクの解法は練習問題としてあった。そこでナンバーリンクの解法を自分で記述し、走行させてみた。けなげに、高速に走ってくれた。

ここでは、その使用経験を紹介することにする。以下の目次に従って紹介する。なお当方、線形計画法や整数計画法はずぶの素人であり、頓珍漢なことを主張しているかもしれない。その点は、ご容赦ください。

 

(1) 線形計画法、整数計画法

(2) GLPKGNU MathProg言語

(3) 制約条件の線形表現

(4) 魔方陣の解法

(5) 数独の解法

(6) ナンバーリンクの解法

 

現在(130926)15個程度のペンシルパズルを整数計画法(mathprog)で攻略している。「キューブ王」または「キューブ王 海永」で当方のトップページ(知恵の輪パズル、キューブパズル 攻略法(解法))googleサーチし、そこからlinkを辿れば、それらの攻略ページに到達できるであろう。そのうち1つが、このページである。

 

 参考にした記事、またはglpkダウンロードのための記事の内容などをいくつか紹介しておく。

記事[]は、GLPK(GNU Linear Programming Kit)を使用し、さらにそれで使えるCPLEX LPという言語のコードをCコードで生成する想定になっていた。しかしここでは別のやり方を採る。GLPKを使用するのは同じ。AMPL(or GNU MathProg)言語で直接に書くものとする。

なおGLPKは無料のツールである。ダウンロードは[]などからたどればいい。スーパー簡易マニュアルというのも用意されている(日本語のマニュアル)。また、ダウンロードしたフォルダ内にdocフォルダがあり、そこに英文のマニュアルがある。コード例はexamplesフォルダ内にある。

なお[2]は、無料ホームページサービスが終了するらしく、消えて無くなるらしい。同じグループの人達が運営すると思われるページを見つけた。それが[3]cygwinなら、そこの「GLPK ダウンロード Full」でいい。[3]も見えなくなったようで、[4]を見つけた。これは、glpkのページなのだろう(英語)。ミラーからのダウンロード、解凍、configure,make、インストールという手順を踏んでいく必要がある(ファイルINSTALLの指示通りでいい)。なお、cygwinmathパッケージ内のglpkは手軽にインストールできるが、出来上がったglpkが極端に遅いので推奨できない。

 

(1)   http://www.is.titech.ac.jp/~okamoto/PDF/2007/puzzleip.pdf

(2)   http://mukun_mmg.at.infoseek.co.jp/mmg/glpk/first.htm

(3)   応用プログラミング演習2006

(4)   http://www.gnu.org/software/glpk/glpk.html

 

1. 線形計画法、整数計画法

 簡単な線形計画法(LP)の問題を提示する。2つの(実数)変数x,yが以下の制約を満たしているとする。

x+4y<=30     (A)

x+y <=15      (B)

x>=0, y>=0

この制約を満たした範囲で、以下の式を最大化せよ。

       F= 30x+45y

こういう問題を解くのは簡単である。(x,y)の範囲は凸領域である。目的関数Fは線形である。目的関数Fは、凸領域の頂点で最大値を取る。今の場合、頂点は4つしかなく、しかも(A),(B)の交点である頂点(x,y)=(10,5)で最大値を取るに決まっている。そこで、以下となる。

        Max F = 30*10+45*5=525

今は、2変数の簡単な例を示した。実際には、1万変数のつまり1万次元の最大化問題を解かなければならない、ということがあるのだそうだ。すると、解くのが大変になり、解けそうもない、ということになる。しかし、シンプレックス法やカーマーカー法などという優れた解法があり、1万変数でも気楽に解けるのだという。こういう話は、数理計画法(or最適化数学)の教科書には必ず載っている(カーマーカー法の精密な話は載らないだろうが)

 さて次に、整数計画法。変化する範囲が整数に制約された変数がある問題の解法を整数計画法という。全変数が整数なら全整数計画法(IP)という。一部が整数なら混合整数計画法(MIP)という。IPMIPは、LPに余分な制約が加わった訳で、LPよりも簡単になる気もする。しかし、実際には100倍以上難しくなるのだという。100倍以上時間が掛かるということらしい。LPなら1秒で解が求まるのだが、MIPなら1分経っても終わらない。1時間経っても終わらない。そういうことのようである。

 だがしかし、落胆することはない。我々は、パズル解法を問題としていた。この場合、整数が変動する範囲は0/1に限定されるようにもできる。全変数が0/1なら0-1整数計画法というらしい。この場合は、解のサーチに特別なやり方が考えられるそうで、かなりの規模の問題も解けるのだという。特別なやり方、分岐カット法などというのがあるらしいが、詳しい話はまったく知らない(ということにする)。言及しない。

 本記事は、0-1整数計画法を活用してパズルを解くのが目標であるといえる。

 

2. GLPKMathProg言語

 線形計画問題を解く商用ツールがいくつかあるようである。AMPLというのが有名らしい。ベル研で開発されたらしい。これらのツールは、勿論、整数計画問題や0-1整数計画問題を解くこともできる。AMPLなどはかなり高価なのであろう。購入45万、保守10万とか。

 しかし、GNUGLPK(GNU Linear Programming Kit)というのを無料配布している。これは、AMPLほど高速ではないのだが、機能的にはAMPLと同等らしい。学生さんなどは、GLPKを使用して、線形計画法や整数計画法を勉強するものらしい。

 このGLPK、複数種類のモデル記述言語を受け付ける。GLPKを起動するとき”-h”を指定して起動すればhelp情報を表示してくれる。記述言語については以下のようなhelp情報がでる。

 

   --glp             read LP/MIP model in GNU LP format

   --mps             read LP/MIP problem in fixed MPS format (default)

   --freemps         read LP/MIP problem in free MPS format

   --cpxlp           read LP/MIP problem in CPLEX LP format

   --math            read LP/MIP model written in GNU MathProg modeling language

   -m filename            read model section and optional data section from filename (the same as --math)

 

 当方は、インターネットの記事に従い、”-m ファイル名GLPKを起動するようにしている。その場合は、GNU MathProg modeling Languageという言語でモデルを記述することになる。

 この言語での記述例を示しておく[2]1円、10円、25円の3種のコイン何枚かで44円を構成せよ。コインの合計枚数を最小化せよ。以上が問題。この記述例が以下。ファイルcoin_m.txtに作るものとする。

 

## ----------- coin  -----------------        ##(A)

set COIN;         ## c1,c10,c25                   ##(B)

param value{COIN};                                    ##(C)

param total;                                           ##(D)

var X{COIN} ,>=0 ,integer;                             ##(E)

 

minimize Number: sum{c in COIN} X[c];                  ##(MIN)

s.t. Total: sum{c in COIN}value[c]*X[c] = total;       ##(SUB)

 

data;    ## ----------- data ------------               ##(Data)

set COIN:= c1 c10 c25;                                  ##(B1)

param value:=                                           ##(C1)

c1      1

c10     10

c25     25;

param total:=44;                                         ##(D1)

end;

 

少し説明しておく。(A)などに現れた”#”は行末コメント(行末までコメント)。

まず(B)で、コイン集合COINを宣言。その構成要素は(B1)c1,c10,c25と定義している。

次に(C)で、コインの属性としての価値valueを宣言。ここの価値は(C1)で、1円、10円、25円と定義している。

次に(D)で、合計金額totalを宣言。(D1)で具体的には44円と定義している。なおparamCの記号定数相当。

次に(E)で、個別コインの使用枚数Xを宣言。”var”としてあり、これを変化させて最適解を求めることになる。

       勿論、非負(>=0)で整数(integer)という制約付きの変数。

次に(MIN)で、目的関数Numberを定義。最小化(minimize)すべき式は “sum{c in COIN} X[c]”。

次に(SUB)で、制約条件Totalを定義。sum{c in COIN}value[c]*X[c]totalに等しいというもの。

       なお、s.tsuch thatかと思ったのだが、subject toだそうである。

次に(Data)。この行以降に、具体的な定義が置かれるというもの。

以降の(B1),(C1),(D1)は既に説明してある。(C1)の構文は説明しておく。C言語の配列への初期値指定と同等である。ただし、Cだと添え字は0,1,2,と変化していくに決まっていて、添え字は省略できる。しかしGNU MathProgでは、添え字はc1,c10,c25などである。添え字とその値の組で指定しないとならないらしい。

最後のend;、これは当然、モデル記述の完了を示す。

 

さて、GLPKを起動(glpsol -m coin_m.txt -o out.txt)すれば、以下のような解が(out.txt)得られる。

 

Problem:    coin_m

Rows:       2

Columns:    3 (3 integer, 0 binary)

Non-zeros:  6

Status:     INTEGER OPTIMAL

Objective:  Number = 8 (MINimum) 1.76 (LP)

 

   No.   Row name        Activity     Lower bound   Upper bound

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

     1 Number              8

     2 Total               44            44             =

 

   No. Column name       Activity     Lower bound   Upper bound

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

     1 X[c1]        *        4             0

     2 X[c10]       *        4             0

     3 X[c25]       *        0             0

 

End of output

 

問題の名前はcoin_m(実際はファイル名で、”.”以降を除けたもの)

問題の規模row,columは、2*3Row,columnは行列の行と列相当。行数は制約条件数(目的関数も含む)、列数は変数の数。

Non-zerosは、制約条件に現れた項数ということだろう。なお、目的関数も制約条件として数える。

求解の状態(status)は、整数の最適解を得ている、というもの。

目標numberの最小値は、整数問題としては8個、実数問題(LP)としては1.76個。

1円4個と10円4個で44円とするのが、整数問題での最適解。計8個のコインを使用。

251.76個で44円とするのが、実数での線形計画問題の最適解。

なお、”*”が付いたものは整数変数で、activity欄が最適解を与える値である。

 

3. 制約条件の線形表現

 整数計画問題や0-1整数計画問題、普通には制約条件は線形制約でなくてもいいはずだろう。だがしかし、GLPKなどでは、線形制約でないとならないのだという[1]

 パズルの制約などは、if(..) then a=1 else a=0などという制約が普通だと思うが、これを線形制約で表現しないとならない。Trueを1、falseを0と表現するとして、変数a,b,cなどがfalse/true(0/1)しか取らない前提で、a^b,aàbなどをどう表現するかを知っておかないとならないのだという。そういうものをいくつか示しておく。

         a&b ß--à a+b=2   a|b  ß--àa+b>=1   aàb ß--à a<=b     a ßà b ß--à a=b

今のは、a,b0/1しか取らなかったので話は簡単であった。

 m1,m2などが1..5の整数をとるとして、m1<m2 à a や m1<m2 ßà aはどうなるか。なおa0/1しか取らない。

      m1<m2 àß--à m2-m1<=5*a      m1<m2ßàß--à m2-m1<=5*a and m1-m2+1<=5*(1-a)

なお、線形制約の方は、”<”は駄目で”<=”でないとならないそうである。

 

4. 魔方陣の解法

 3*3の魔方陣の例。[1]の例からもってきた。3*3のマス目に1 から9 までの数字を1つずつ入れて、どの横行、縦列の和も同じになるようにするのが問題。

 

ban

11

12

13

 

 

1

5

9

 

21

22

23

 

 

8

3

4

 

31

32

33

 

 

6

7

2

 

3*3の盤面に11,12などの番号を与える(左図)。個別のマスに右図のような数字を与えれば、それが1つの解となる。

この問題の解法は以下のコードでいい。

 

## --------- magic -----------------

set I123:=1..3;

set Dig:=1..9;

var ban{I123,I123,Dig} ,binary;

var Ban{I123,I123}, integer, >=0, <=9;

 

s.t. uniq0{r in I123, c in I123}: sum{n in Dig} ban[r,c,n] =1;

s.t. uniq1{n in Dig}: sum{ r in I123, c in I123} ban[r,c,n] =1;

s.t. row{r in I123}: sum{ c in I123,n in Dig} n*ban[r,c,n] =15;

s.t. col{c in I123}: sum{ r in I123,n in Dig} n*ban[r,c,n] =15;

s.t. conv{r in I123,c in I123}: sum{d in Dig}d*ban[r,c,d] =Ban[r,c];

 

minimize value: ban[1,1,1];

solve;

for{r in I123}{for{c in I123}printf "%2d",Ban[r,c]; printf("\n");}

 

end;

 

banは本来は2次元なのだが、binary(0/1)変数にするために3次元となっている。

制約uniq0は、ban[i,j,*]を1価関数とするもの。制約uniq1は、個別のマスに1..9を1つずつ割り振るもの。

制約rowは、個別の行で3つの値を加えれば15になるというもの。

制約colは、個別の列で3つの値を加えれば15になるというもの。

制約convは、3変数binaryban[,,]を2変数のBan[,]に変換するもの。solve以降のfor()でのダンプ対策。

最小化(minimize)の行は本当はいらない。あれば最適解を探す。なければ、可能解を探して終わる。ここは1つの可能解でいい。

次のsolve,for{}、この2つ、マニュアルをみた範囲では、そういう記述はなかった。しかし、ダウンロードしたGLPKexamplesフォルダ内のthai.modファイル内などにあった。solveは当然、「解け」ということだろう。これがあると、以降のfor}が撥ねつけられなくなる。このfor{}は勿論、3*3の盤面状態を3*3の形にダンプするもの。

 

結果(out.txt)を抜粋すると以下。ban[r,c,d]の内、trueとなったもののダンプ。

 

ban[1,1,2]   *              1             0             1

 ban[1,2,6]   *              1             0             1

 ban[1,3,7]   *              1             0             1

 ban[2,1,4]   *              1             0             1

 ban[2,2,8]   *              1             0             1

 ban[2,3,3]   *              1             0             1

 ban[3,1,9]   *              1             0             1

 ban[3,2,1]   *              1             0             1

 ban[3,3,5]   *              1             0             1

 

説明不要であろう。説明はしない。

上のコードの太字部分でのダンプ結果は以下となる。標準出力に吐き出される。

 

2 6 7

4 8 3

9 1 5

 

 

5. 数独の解法

 魔方陣を少し難しくしたのが数独。魔方陣でなく、ラテン方陣を少し難しくしたものというのが正しいか。

 Wikipediaでのラテン方陣の説明を載せておく。

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

ラテン方格(-ほうかく、Latin square)とは n n 列の表に n 個の異なる記号を、各記号が各行および各列に1回だけ現れるように並べたものである。ラテン方陣(-ほうじん)ともいう。例を示す:

 

1

2

3

2

3

1

3

1

2

 

ラテン方格は数学的には擬群の積表と見ることができる。

ラテン方格は実験計画法に応用される。またペンシルパズルの一種「数独」もラテン方格の応用である。

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

数独は9*9が標準。使う数字は1..920-30個のマスには予め数字が入っている。残りのマスに数字を入れてラテン方陣を完成する。さらに9個の3*3小方陣にも1..9の数字を1つずつ入れる制約がある。

 

0

0

5

3

0

0

0

0

0

 

1

4

5

3

2

7

6

9

8

8

0

0

0

0

0

0

2

0

 

8

3

9

6

5

4

1

2

7

0

7

0

0

1

0

5

0

0

 

6

7

2

9

1

8

5

4

3

4

0

0

0

0

5

3

0

0

 

4

9

6

1

8

5

3

7

2

0

1

0

0

7

0

0

0

6

 

2

1

8

4

7

3

9

5

6

0

0

3

2

0

0

0

8

0

 

7

5

3

2

9

6

4

8

1

0

6

0

5

0

0

0

0

9

 

3

6

7

5

4

2

8

1

9

0

0

4

0

0

0

0

3

0

 

9

8

4

7

6

1

2

3

5

0

0

0

0

0

9

7

0

0

 

5

2

1

8

3

9

7

6

4

 

左が初期状態(0は実際は未定)。右が答え。最初23個しか数字が確定してないので、かなり難しい問題である。

コードは以下。

 

## ----------- sudoku -------------------

set I19 := 1..9;        # row/column index

set D19 := 1..9;       # digits

set I13 := 1..3;        # mat index

set I036:= 0..6 by 3;              # mat offset

param ini{I19,I19} ,integer;      # pre-located digits

var ban{I19,I19,D19} ,binary;   # board[r,j] is d

var van{I19,I19} ,integer, >=1, <=9;

 

s.t. uniq{ r in I19, c in I19}: sum{ n in D19} ban[r,c,n] =1;       # ban[r,c]=n

s.t. hori{ r in I19, n in D19}: sum{ c in I19} ban[r,c,n] =1;        # Latin square

s.t. vert{ c in I19, n in D19}: sum{ r in I19} ban[r,c,n] =1;        # Latin square

s.t. mat { r in I036, c in I036, n in D19}:                              # Sudoku square

                             sum{ i in I13, j in I13} ban[r+i,c+j,n]   =1;

s.t. init{ r in I19, c in I19,       n in D19: n=ini[r,c]}:               ban[r,c,n] =1;

s.t. conv{ r in I19, c in I19}: van[r,c]=sum{ n in D19} n*ban[r,c,n];

 

minimize value: ban[1,5,1]+5*ban[9,9,1];

solve;

for{r in I19}{ 

#             for{c in I19} printf " %s", sum{n in D19}n*ban[r,c,n];

              for{c in I19} printf " %s", van[r,c];

              printf("\n");

}

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

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

              1             0 0 5 3 0 0 0 0 0

              2             8 0 0 0 0 0 0 2 0

              3             0 7 0 0 1 0 5 0 0

              4            4 0 0 0 0 5 3 0 0

              5             0 1 0 0 7 0 0 0 6

              6             0 0 3 2 0 0 0 8 0

              7             0 6 0 5 0 0 0 0 9

              8             0 0 4 0 0 0 0 3 0

              9             0 0 0 0 0 9 7 0 0;

#            1             0 2 0 0 7 5 0 1 0

#            2             1 0 0 0 0 4 5 0 8

#            3             0 5 6 8 0 0 2 0 0

#            4            8 1 0 0 0 0 7 0 0

#            5             9 0 0 0 0 0 0 0 3

#            6             0 0 2 0 0 0 0 4 5

#            7             0 0 9 0 0 1 4 3 0

#            8             3 0 1 7 0 0 0 0 9

#            9             0 7 0 3 9 0 0 2 0;

end;

 

特別に説明するほどのことはない。ただし2つ指摘しておく。

最小化(minimize)の指定はなくていい。あると、最適解を探して遅くなる。ないと、可能解を1つ見つけると終わる。結局速くなると書いてある記事を見た覚えがある。しかし、速くなるとは限らないようである。

 走らせると、正しい答えを見つけてくれた。当方の普通のPCだと、0.6秒程度で終わった。

 

6. ナンバーリンクの解法

 当方、実は数独はあまり評価してない。簡単な問題なら、手間は掛かるかもしれないが必ず正解に到達する。近頃、バックトラックを必要とする難しい問題が流行っているようだが、これだと正解に到達しそうもない。いずれにしても、やり終わって達成感というか感動を覚えることはない。

 一方、ナンバーリンクは好きである。大昔、二コリのペンシルパズル本シリーズを眺めていて、ナンバーリンクが面白そうなので買った。家に帰ってやってみると、かなり難しかった。悪戦苦闘して正解を見つけ、正解の線が大きく迂回している部分などを観賞し、感動を覚えていた。直ぐにナンバーリンクのファンになったのだが、このナンバーリンクその当時、1冊しか販売されてなかった。残念に思った。今の時点では3冊目まで出ている。数独は数え切れない程出ているのに、肝心のナンバーリンクが3冊である。早く4冊目が出て欲しい。

 数独は、計算機で解くのも簡単だし、問題を作るのも簡単らしい。一方、ナンバーリンクは、解くのは難しい。人が解くのも難しいし、計算機で解くのも難しいらしい。この辺りは、後で触れる。

 

 資料[1]では、ナンバーリンクのモデルの捕らえ方についての説明がない。それで、自分で考えた。それを示しておく。

 wikipediaの記事にあった3*4行列で説明する。個別の要素位置をマスと呼ぶ。初期状態では、いくつかのマスに数字が割り振られる。1が2つ、2が2つなどである()。最初から数字が割り振られたマス位置を終端マスと呼ぶ。そうでないマスを途中マスと呼ぶ。課題は、同じ数字を線で結ぶこと(交差は不可)。中の2つの図がその例。課題を作る立場では、複数の結線群(つまり複数の解)が存在してはならない。結線群から置いてきぼりされたマス(右から2番目)があってはならない。なお、2番目の制約は、ニコリでは採用されている。2番目の制約を採用しない一派もいる。置いてきぼりされたマスがある解を短絡解または関西解というが、こういう解を最初に指摘したのが関西の人で、この解を気にいっているのが関西の人だという。関西解を許す課題は、極端に難しくなるのだそうだ。また、関西解がないことを確かめるのも極端に難しいという。最初の制約、つまり解が1つという制約は、関西流でも採用されているようだ。

さて、人間が鉛筆で解を探していく場合は、結線していくのがやり易い。しかし、計算機で解くとなると別の考えがいいようである。

右の図を考える。個別のマスの左右上下にドアがあるとする。ドアは開いている(|,-)か閉じている()。終端マスは1つのドアが開いている。途中マスは、置いてきぼりされてないのなら、2つのドアが開いている。置いてきぼりされたマスは4つのドアが全て閉まっている。

ドアが開いて2つのマスがつながっていると、同じ数字が割り当てられる。また、角のマスは外に接する2つのドアは閉まっている。辺のマスは外に接する1つのドアは閉まっている。右上の角マスは、終端マスではない。今は4つのドアが全て閉まっている。しかし、関西解を許さないなら、2つのドアは開いているはずである。左下の角マスも同様。

中の2つのマス、一方は2で他方は1である。共有のドアが開いていてはならない。開いていれば同じ数字でなければならないからである。

 

さて、以上の考えに従う7*7の盤面での求解コードを示す。8*8にするには、Nr,Ncを変更要。さらに、使う数字が6以上になればNdも変更要。

 

## -------------------- number link ------------------------

param Nr:=7;                                       ## number of rows

param Nc:=7;                                       ## number of columns

param Nd:=6;                                      ## number of degits

set Ir:= 0..Nr;                                      ## row indexes

set Irn:= 0..Nr-1;                  ## row indexes, narrow

set Ic:= 0..Nc;                                     ## column indexes

set Icn:=0..Nc-1;                   ## column indexes, narrow

set In :=1..Nd;                                     ## number

param ini{Irn,Icn},integer;        ## initial board

var hs{Irn,Ic} ,binary;                            ## adjacent room has same digit(horizontal)  (H

var vs{Ir,Icn} ,binary;                             ## adjacent room has same digit(vertical)   (V)

var b{Irn,Icn,In} binary;                         ## ban[r,c] == n                               (B)

 

#minimize van: b[2,2,1];          ## objective [dummy] function

 

s.t. unique{ r in Irn, c in Icn}: sum{ n in In} b[r,c,n] =1;           # ban[r,c]=n

s.t. vvv0{ c in Icn}: vs[0,c]=0;   # wall

s.t. vvv1{ c in Icn}: vs[Nr,c]=0; # wall

s.t. hhh0{ r in Irn}: hs[r,0]=0;   # wall

s.t. hhh1{ r in Irn}: hs[r,Nc]=0; # wall

s.t. init1{ r in Irn, c in Icn, n in In: n =ini[r,c]}: b[r,c,n]=1;

s.t. same0{r in Irn,c in Icn diff{0},n in In}:

              hs[r,c]      <=b[r,c-1,n]-b[r,c,n]+1;   #hs[r,c]=1-->b[r,c,n]=b[r,c+1,n];

s.t. same1{r in Irn,c in Icn diff{0},n in In}:

              hs[r,c]      <=b[r,c,n]-b[r,c-1,n]+1;   #hs[r,c]=1-->b[r,c,n]=b[r,c+1,n];

s.t. same2{r in Irn diff{0},c in Icn,n in In}:

              vs[r,c] <= b[r-1,c,n]-b[r,c,n]+1;  #hv[r,c]=1-->b[r,c,n]=b[r+1,c,n];

s.t. same3{r in Irn diff{0},c in Icn,n in In}:

              vs[r,c] <= b[r,c,n]-b[r-1,c,n]+1;  #hv[r,c]=1-->b[r,c,n]=b[r+1,c,n];

s.t. open1{r in Irn,c in Icn:ini[r,c]>=1}: # terminal rooms have one open door

                             vs[r,c]+hs[r,c]+vs[r+1,c]+hs[r,c+1]=1;

s.t. open2{r in Irn,c in Icn:ini[r,c]=0}:  # mid rooms have two open doors

                             vs[r,c]+hs[r,c]+vs[r+1,c]+hs[r,c+1]=2;

## -------------------- data ------------------------

data;

param      ini:          0 1 2 3 4 5 6 :=

              0             1 0 2 3 0 0 1

              1             0 0 0 0 0 3 0

              2             0 2 0 4 0 0 0

              3             0 0 0 0 0 0 0

              4             4 5 0 0 0 5 0

              5             0 0 0 0 0 0 0

              6             6 0 0 0 0 0 6          ;

end;

 

最初に、変数の説明。(H),(V),(B)の行。

まず(B)。b[ , , ]7*7*6である。7行で7列、可能な数字が1..6ということ。

次に(H)hs[ , ]7*8である。行は7、列位置のドア数は植木算で7+=8となる。1はsame or openを意味し、0don’t care or closedを意味する。

次に(V)。これは、行と列の立場が逆転しているが、(H)と同様。省略。

マスに対応するのがb[ , , ]。行と列の添え字とそこに割り振られた数字を表現する。

 さて次に、制約条件の説明。

制約uniqueは、各マスに割り当てられた数字は1つという制約。

制約vvv0,vvv1,hhh0,hhh1は、外部に接したドアは閉まっているという制約。

制約init1は、終端マスの設定。

制約link1は、終端マスなら1つのドアが開いているという制約。

制約link2は、途中マスなら2つのドアが開いているという制約。結局、関西解は許してないことになる。

制約same0,sam1は、横方向のドアが開いていれば、つながった2つのマスは同じ数字となるという制約。

制約same2,sam3は、縦方向のドアが開いていれば、つながった2つのマスは同じ数字となるという制約。

 

以上のコードでも、7*78*8程度のナンバーリンクは求解できた。1秒以内である。それではと10*10の問題をやってみた10分以上たっても終わらなかった。この点は、予め類推できていたともいえる。

 実は、ナンバーリンクの課題、関西解は許さないという約束なら、解く場合には極端に強力な「角の定理」が使用できる。

 

 

結線で直角が確定するなら、直角が定める斜め方向に直角が伝播する。終端マスにぶつかるまで伝播する。というもの。たとえば、マスa1は直角の結線で確定である。すると、a2,a3,a4,a5,a6も直角の結線となる。またb1で直角の結線が確定したとする(7行程度下の太字部分参照)。するとb2,b3,b4,b5も直角の結線である。

 個別のドアの開閉を正しく指定できたとすれば、正しく結線できるだろう。ドアの開閉未定のドア数は7*7の問題なら、30個程度か。つまり2^30程度の場合を調べればいい。しかし10*10の問題なら、ドアの開閉未定のドア数は80個程度か。なら2^80程度の場合を調べなくてはならない。終わる筈がない。

 しかし、角の定理が使用できれば、場合の数を極端に削減できる可能性がある。2^802^30になれば終わってくれる期待が持てることになる。角の定理の他に「辺の定理」というのもありそうで、そうであれば、10*10問題も気楽に解ける可能性がある。個人的には10*18程度は気楽に解いてくれないとならない、と思っている。

 なお、この問題では数字1は、辺を這っていく必要がある。辺の定理の一種。するとb1の直角が確定する。

 角の定理や辺の定理を導入するつもりにはなっているのだが、まだ導入してない。代わりに、「Uターン禁止の定理」を導入してみた。個別のマスの4隅に柱があるとして、柱を共有する4つのマスを考える。柱の傍の4つのドアの内3つが開いていれば、結線がUターンする。それで解があるとする。しかし、残りの1つのドアを開ければUターン部分を通らないで、関西解での結線が可能となる。そこでUターンは許さない。それが「Uターン禁止の定理」。

 

 

この効果を示しておこう。

 

 

 左の図、外に面したドアを全て閉じている。また、関西解を拒否することで右上と左下のマスの内側の赤いドアが開いている(赤い1)。次に左から2番目の図を考える。Uターン禁止から青色のドアが閉まることになる(青い丸)。次に右から2番目の図を考える。右や左の辺のマスは途中マスでありピンクのドアが開くことになる(ピンクの1)。以降も同様な論理展開で、右端の結線となろう。つまり、解が得られたことになる。

 これに従うコードが以下。盤面は3次元のb[ , , ]から2次元のban[ , ]に変更した。この点注意。

 

## -------------------- number link ------------------------

param Nr:=10;                                                   # number of rows

param Nc:=18;                                                   # number of columns

param Nd:=9;                                                     # number of digits

#param Nr:=14;                                                  # number of rows

#param Nc:=14;                                                 # number of columns

#param Nd:=13;                                                 # number of digits

set DIG:=1..Nd;                                                 # 1 2 3 4 5 6 7 8 9 a b ...

param ini{0..Nr-1,0..Nc-1},integer;         # initial board

param aka{0..Nr-1,0..Nc-1},default 0,binary;          # 0/1: normal/hidden

var hd{0..Nr-1,0..Nc},binary;                 # 0/1: closed/open (horizontal)

var vd{0..Nr,0..Nc-1},binary;                 # 0/1: closed/open (vertical)

var ban{0..Nr-1,0..Nc-1},integer,>=1,<=Nd;           # ban[r,c]

 

s.t. vvv0{ c in 0..Nc-1}: vd[0,c]=0;          # outer doors are closed

s.t. vvv1{ c in 0..Nc-1}: vd[Nr,c]=0;        # outer doors are closed

s.t. hhh0{ r in 0..Nr-1}: hd[r,0]=0;          # outer doors are closed

s.t. hhh1{ r in 0..Nr-1}: hd[r,Nc]=0;        # outer doors are closed

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

              hd[r-1,c]+hd[r,c]+vd[r,c-1]+vd[r,c]<=2;  # inihibit U-turn

s.t. open2{r in 0..Nr-1,c in 0..Nc-1:ini[r,c]=0}:

#             vd[r,c]+hd[r,c]+vd[r+1,c]+hd[r,c+1]+2*aka[r,c]=2;#mid room: 2 open doors

              vd[r,c]+hd[r,c]+vd[r+1,c]+hd[r,c+1]=2;    #mid room: 2 open doors

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

              vd[r,c]+hd[r,c]+vd[r+1,c]+hd[r,c+1]=1;    # terminal room: 1 open door

s.t. close_h{r in 0..Nr-1,c in 0..Nc-2:ini[r,c]>=1 and ini[r,c+1]>=1}:

              hd[r,c+1]=0;                                        # ini[r,c]!=ini[r,c+1] --> hd[r,c+1]=0

s.t. close_v{c in 0..Nc-1,r in 0..Nr-2:ini[r,c]>=1 and ini[r+1,c]>=1}:

              vd[r+1,c]=0;                                        # ini[r,c]!=ini[r+1,c] --> vd[r+1,c]=0

 

s.t. init1{r in 0..Nr-1,c in 0..Nc-1,n in DIG:n =ini[r,c]}: ban[r,c]=n;#terminal

s.t. same0{r in 0..Nr-1,c in 1..Nc-1}: ban[r,c]-ban[r,c-1]<=(Nd-1)*(1-hd[r,c]);

s.t. same1{r in 0..Nr-1,c in 1..Nc-1}: ban[r,c-1]-ban[r,c]<=(Nd-1)*(1-hd[r,c]);

s.t. same2{r in 1..Nr-1,c in 0..Nc-1}: ban[r,c]-ban[r-1,c]<=(Nd-1)*(1-vd[r,c]);

s.t. same3{r in 1..Nr-1,c in 0..Nc-1}: ban[r-1,c]-ban[r,c]<=(Nd-1)*(1-vd[r,c]);

 

#minimize objFunc: 3*ban[2,2]+ban[4,4];               # objective dummy function

solve;

for{r in 0..Nr-1}{ for{c in 0..Nc-1} printf "%3d", ban[r,c]; printf("\n");}

data; ## --------------- problem -------------------

#param aka[1,0]:=1;

param      ini:          0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17:=

              0             1  0  0  0  0  0  0  0  0  0  0  0  0  0  3  0  0  0

              1             0  0  2  0  0  0  0  0  0  0  0  4  0  0  0  0  9  0

              2             0  0  0  0  0  0  7  0  0  0  0  0  0  6  0  0  0  0

              3             0  0  3  0  0  0  0  0  0  0  0  0  0  0  0  0  8  0

              4             0  0  0  0  0  0  0  0  0  0  0  0  7  0  0  0  0  0

              5             0  0  0  0  0  0  0  0  6  0  0  0  0  0  0  0  0  0

              6             0  0  0  0  0  0  0  0  0  0  0  0  5  0  0  0  0  0

              7             0  0  4  0  0  0  0  0  0  0  0  0  0  0  9  0  0  0

              8             0  0  0  0  5  0  8  0  0  0  0  0  0  0  0  2  0  0

              9       0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1;

#param    ini:          0  1  2  3  4  5  6  7  8  9 10 11 12 13:=

#            0             1  0  0  0  0  0  0  0  0  0  0  0  0 11

#            1             0  0  3  0  0  4  5  6  0  0  0  0  0  0

#            2             0  0  0  0  0  0  0  0  7  5  0  0  0  0

#            3             0  0  8  0  0  0  0  0  0  0  0  0  0  0

#            4             0  0  9  0  0  0  0  0  0  0  0  0  0  0

#            5             0  0  0 10  0  0  0  0  0  6  0  0  0  0

#            6             0  0  0  0  0  0  0  0  0  0  0  0  0  0

#            7             0  0  0  0  0  5  0  0  7  0  0  0  0  0

#            8             0  0  0  0  0  0  0  0  0  0  0  0  0  0

#            9      13  0  0  8  0  0  0  0 11 12  0  0  0  0

#            10           0  9  0  0  0  0  0  0  0  0  0  0  0  0

#            11           0  0  0  1  0  0  0  0  0  0  0  0 10  0

#            12           0 13  0  0  0  0  0  5  0  0  0  0  0  0

#            13           2  0  0  0  0  0  0  2  3  4  0  0  0 12;

end;

 

制約条件は、いくつか指摘しておく。

制約notUは、Uターン禁止の制約。

制約same0,same1は、横方向のドアが開いていれば、つながった2つのマスの数字は同じという制約。

制約same2,same3は、縦方向のドアが開いていれば、つながった2つのマスの数字は同じという制約。

14*14問題を解けた。1.6秒。10*18も解けた。1.7秒。

 

このコードは、glpkの版によって、実行時間がかなり変わる。0章の「導入」部で説明した[2],[3]のものであれば、1.6秒程度である。しかし[4]であると、10秒以上かかった(変数banintegerを指定した場合)。しかし、banintegerを除ければ(banを実数変数とすれば)0.5秒程度となり、速くなる。一般には、[4]の方が[2],[3]よりも数倍速いようである。

 細かい話をしておく。 [4]を作る(make)際、"make"だけであると、gccの選択子が"-O2 -g"となるようで、"make CFLAGS=-O3"とすればgccの選択子が"-O3"となるようである。実行速度は大して変わらないようだが、ロードモジュールのサイズが半分になる。

 

関西解をサポートする手当ても一部入れようとした。Param aka[ , ]は開かずの間または隠された部屋を意味する。これが1なら、そのマスのドアは総て閉まっている。

 なお、途中マスが総て開かずの間になる可能性があるとして、つまりaka[ , ]を変数として、しかも0/1どちらも可能として、解をサーチさせると当然極端に遅くなる。7*710分程度は走行させたことがある。終わらなかった。

 開かずの間となる可能性があるマスを上手に特定し、それで走らせると、解を求める可能性はあろう。しかし、今の時点では、開かずの間となる可能性についての知識は何もない。結局、関西解はサポートしてない。

 

トップページへ

 

Ads by TOK2