GLPKでのスリザーリンクの攻略法(解法)

キューブ王

最初のコード(111018)は、速度的に少し不満であった。そこで今回コードを調整した(130415)。記載内容も少し修正した。

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

パズルの問題を解くプログラムについて挑戦している。今回、スリザーリンクに挑戦。GLPKMathProgで書いた。

このスリザーリンク、「1つのリングになること」という制約条件の指定が難しい。実は、独力では指定法を編み出せなかった。しかし、整数計画法で有名な巡回セールスマン問題の解法を活用できることが分かった。そして、実際にMathProgで書いてみたもの。少し遅いかとは思うが、10*10問題なら許容できるであろう。

 

[スリザーリンクのルール]      (図はインターネットでサーチした杉村氏の記事を参考にした)

 

 

一般にはNr*Ncだが、4*4の例でしめす。4*4個のマスで構成される盤面がある。いくつかのマスには0-3の数字が書いてある。

 課題は、マスの辺に沿って線を引いていくもの。

(1)   1つのリングを作る。2つなど複数のリングになると駄目。

(2)   数字が書いてあるマスはその数字と同じ数の辺を通過する。

(3)   頂点を共有する辺は4個(隅などでは2-3となる)であるが、通過した辺は0or2であること。1,3,4個は駄目。

上の図、リングの外を黄色、リングの中を灰色で色分けしてある。

左の盤面、正解となる。真ん中の盤面は、リングが2つあり駄目。右の盤面は、通過した辺が4,3,1の頂点があり駄目。

 なお、「スリザー」とはslitherで、英和辞典で調べると、「ずるずる[するする]滑る、〔蛇のように〕滑るように進む」とある。結線を蛇の胴体や尻尾と思えば、1つのリングとなっている蛇であり、ウロボロスということになる。また、ある頂点で3つに別れた結線を考えると、その内の1つは足であろうから、それを禁止することは「蛇足」を禁止することともいえよう。

 なおリングの中を灰色で表示する前提で、さらに蛇足が現れない前提では、答えの妥当性を確認するための表示で(太い黒で)線を引く部分を省略できる。灰色領域と黄色の領域との境界線が、丁度太い黒線部分に対応するからである。以降、「太い黒線」を省略することもあろう。

 

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

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

 

 

本来の盤面は、黄色いマス4つ(2*2)である。ノード(小さい四角)は植木算で9(3*3)である。横の辺(横長の四角)や縦の辺(縦長の四角)は左右や上下に伸ばし、12(3*4,4*3)である。なお、マスは右下方(肌色部分)に1つ伸ばし9個(3*3)と考える場合もある。個別の片に付けた1,1等の数値の組みは、個別片を識別する座標である。

 次に本論の制約条件。これは1つを除いて非常に簡単。1つとは、「複数のリングは駄目」のことだが、これは後で説明する。

 

 

まず、数字が指定されたマスでの制約。数字nが指定されていて、マスの辺がa,b,c,dとする(左図)。さらに通過する/しないを1/0で表現しているとする。以下でいい。左は抽象論。右は具体論。

          a+b+c+d=n                 hl[r,c]+vl[r,c]+hl[r+1,c]+vl[r,c+1] = ini[r,c]

次に、リングを形成するという制約。個別の頂点について、その頂点をを共有する辺をa,b,e,fとする。またその頂点を通過してないかどうかをjで表すとする。以下でいい。

         a+b+e+f  +2*j  =2                hl[r,c]+vl[r,c]+hl[r,c-1]+vl[r-1,c]  +2*j[r,c]  =2

これでいいのだが、j[r,c]は変数となる。すると遅くなる。別のやり方がある。これはインターネットで見つけた(杉村氏)。以下でいい。左の方(抽象論)だけ示す。

         a+b+e+f <= 2                         通過は2つまで

         a<=b+e+f, b<=a+e+f, e<=a+b+f, f<=a+b+e  1つ通過していれば後1つ通過。1つだけの通過を拒否

ある頂点を通過しているとすれば、入るのと出るのと2つあるわけで、結局、通過する辺の集まりはリングを構成することになる。

 正式のコード例も示しておく。

 

## -------- must be rings ----------------------------

        var hl{1..Nr+1,0..Nc+1} binary;          # 0/1: not/link (---)

        var vl{0..Nr+1,1..Nc+1} binary;          # 0/1: not/link ( | )

s.t. RNG0{ r in 1..Nr+1}: hl[r,0]=0;     # outer segments

s.t. RNG1{ r in 1..Nr+1}: hl[r,Nc+1]=0;  # outer segments

s.t. RNG2{ c in 1..Nc+1}: vl[0,c]=0;     # outer segments

s.t. RNG3{ c in 1..Nc+1}: vl[Nr+1,c]=0;  # outer segments

s.t. RNG4{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c-1]+vl[r-1,c]+vl[r,c]+hl[r,c]<=2;  

s.t. RNG5{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c-1]<=vl[r-1,c]+vl[r,c]+hl[r,c];    

s.t. RNG6{ r in 1..Nr+1, c in 1..Nc+1}:vl[r-1,c]<=hl[r,c-1]+vl[r,c]+hl[r,c];    

s.t. RNG7{ r in 1..Nr+1, c in 1..Nc+1}:vl[r,c]<=hl[r,c-1]+vl[r-1,c]+hl[r,c];    

s.t. RNG8{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c]<=hl[r,c-1]+vl[r-1,c]+vl[r,c];    

s.t. LNK3{r in 1..Nr,c in 1..Nc:ini[r,c]<=3}:    # 0,1,2,3: number of links

        hl[r,c]+vl[r,c]+hl[r+1,c]+vl[r,c+1]=ini[r,c];

s.t. LNK8{r in 1..Nr,c in 1..Nc:ini[r,c]=8}:     # 8: don't care

        hl[r,c]+vl[r,c]+hl[r+1,c]+vl[r,c+1]<=3;

s.t. LNK9{r in 1..Nr,c in 1..Nc:ini[r,c]=9}:     # 9: "3" with starting node

        hl[r,c]+vl[r,c]+hl[r+1,c]+vl[r,c+1]=3;

 

制約RING0-3は、本来の盤面よりも外の辺はリンク(結線)してない、ということを示すもの。

制約RNG4は、個別ノードで考えて、結線数は2以下というもの。

制約RNG5-8は、個別ノードで考えて、1つ結線していれば、後1つ結線している、ということを示す。

制約LNK3は、問題定義配列で指定された数字ini[r,c]が3以下なら、対応するマスの4辺の結線数もini[r,c]になることを示す。

制約LNK8は、問題定義配列で指定された数字ini[r,c]が8なら、対応するマスの4辺の結線は3以下ということを示す。”8”” “の代用。

制約LNK9は、ここでは”3”の指定と同じ制約となる。別のところでは、”9”は特別な意味を持つ(こともある)

 

[1つのリングとなる制約の記述法1] 小さいリングの拒否

 最後に、1つのリングとなる制約。この制約を指定するのは難しい。Cで書くのなら問題ないはずだが、整数計画法の制約条件として書くのは難しい。独力でかなり考えたが、いい手は編み出せなかった。それではと、インターネットでサーチした。いい手はサーチできなかった。ただ、以下のやり方(切除平面法というらしい)を見つけた(杉村氏)

小さいリングを構成する辺の集まりをピックアップして、その小さいリングを拒否する。具体的には、構成辺を加算する。その加算結果が辺の数よりも1つは少ないという制約を追加するもの。上の図ではa+b+c+d+e+f<=5でいい。 

このやり方で、確かに上手くいく。ただし、(我々の場合)手作業が発生する。そこで、猿知恵ではあるのだが、別のやり方も考えた。1つ目は以下。インターネットの記事のやり方を参考にしたもの。

 

## -------- inhibit small ring, directly -------------------

s.t. SR0{ r in 1..Nr, c in 1..Nc-1}:     #two cells which are holzontary adjacent

        vl[r,c]+vl[r,c+2]+hl[r,c]+hl[r,c+1]+hl[r+1,c]+hl[r+1,c+1]<=5;

s.t. SR1{ r in 1..Nr-1, c in 1..Nc}:     #two cells which are vartically adjacent

        hl[r,c]+hl[r+2,c]+vl[r,c]+vl[r+1,c]+vl[r,c+1]+vl[r+1,c+1]<=5;

 

 隣接した2つのマスだけでリングとなっている小さいリングを最初から明示的に拒否するもの。制約SR0は横長の長方形、制約SR1は縦長の長方形を拒否するもの。なお、1つのマスだけでリングを作るものは以下の制約で初めから拒否してあった。

 

s.t. LNK8{r in 1..Nr,c in 1..Nc:ini[r,c]=8}:     # 8: don't care

        hl[r,c]+vl[r,c]+hl[r+1,c]+vl[r,c+1]<=3;

 

このやり方、少しは効果があった。ただし、マス数3以上のリングを拒否しないので、これで総て解決ということにはならない。

 

[1つのリングとなる制約の記述法2] 目標を定める

2つ目。整数計画問題の目標を定める。通過した縦の辺数と通過した横の辺数を最大化する目標とした。

 

s.t. Length: length= sum{r in 0..Nr+1,c in 0..Nc+1}(hl[r,c]+vl[r,c]);

maximize maxmax: length;

 

少しは、上手くいった気はする。しかし、最終的には却下した。そして、別のやり方を考えた。勿論これも猿知恵でしかない。

 平面上のリングは、領域を2つに分ける。リング内とリング外の2つの領域に分ける訳である。リング外を白(0)、リング内を黒(1)と塗り分けるとする(リング内のリングは白とする)。例えば以下のよう。

 

 

左は、制約条件を満たしてないもの。左上に小さい島があり、また左の中程に池もある。島の周囲を巡る黒線と池の周囲を巡る黒線が小さいリングである。右は、制約条件を満たしているもの。島も池もなく、確かに黒線が、1つのリングになっている。

 島や池を、明確に拒否できればいいのだがそうもいかない。しかし、黒線の数を多くする方がいい気がする。同様に、黒いマスの数を多くする方がいい気もする。この節の最初のやり方、黒線だけを考えたのだが、黒マスも考慮にいれるというもの。

 そういうことで、黒線の数と黒いマスの数を数えることにした。

 

まず、マスを色分けする部分が必要になる。

 

## ---------------- color --------------------------------

        var col{1..Nr,0..Nc} >=0,<=1;    # 0/1: white/black

s.t. COL0{ r in 1..Nr}: col[r,1]=vl[r,1];       

s.t. COL1{ r in 1..Nr, c in  2..Nc}: col[r,c]-col[r,c-1]<=vl[r,c]; # vl=0-->same

s.t. COL2{ r in 1..Nr, c in  2..Nc}: col[r,c-1]-col[r,c]<=vl[r,c]; # vl=0-->same

s.t. COL3{ r in 1..Nr, c in  2..Nc}: vl[r,c]<=col[r,c-1]+col[r,c]; # vl=1-->not

 

制約COL0は、左端のマスの色を設定。左端の縦の線vl[,]が黒/白なら、マスの色col[r,c]を黒/白とするもの。

制約COL1,COL2は、横方向に隣接する2つのマスの境界線vl[,]が白なら、2つのマスの色が等しいとするもの。

制約COL3は、横方向に隣接する2つのマスの境界線vl[,]が黒なら、2つのマスの色が違うとするもの。

 

 そして、黒いマスの数と黒い線の数を数えた。

 

        var area;       # number of cells with col=1

        var length;     # number of links

s.t. AREA:   area  = sum{r in 1..Nr,c in 1..Nc}col[r,c];

s.t. LENGTH: length= sum{r in 1..Nr+1,c in 1..Nc+1}(hl[r,c]+vl[r,c]);

 

 最後に、最大化する式。いくつか用意した。

 

maximize maxmax: length+0.01*area;       # (M1)

#maximize maxmax: length-0.01*area;      # (M2)

#maximize maxmax: -length;               # (M3)

#maximize maxmax: area+0.01*length;      # (M4)

#maximize maxmax: area-0.01*length;      # (M5)

#maximize maxmax: -area;                 # (M6)

 

(M1)では、lengthを最大化し、それからareaを最大化することになる。0.01*areaはほぼ[0,1]に収まる点に注意。

(M2)では、lengthを最大化し、それからareaを最小化することになる。

(M3)では、lengthを最小化することになる。

(M4)では、areaを最大化し、それからlengthを最大化することになる。0.01*length[0,1]に収まる点に注意。

(M5)では、areaを最大化し、それからlengthを最小化することになる。

(M6)では、areaを最小化することになる。

 

 そして、それらの効果を評価してみた。どれか1つで全て旨くいく、ということにはならなかった。1つでやって旨くない時、別のをやれば旨くいく。そういうもの。4種の目標(M1,M2,M4,M5?)を試行すれば、10個程度の問題全てを正解できた()

 

[巡回セールスマン問題]

 巡回セールスマン問題というのがある。複数の町を1人のセールスマンが巡回していく。個別の町に必ず立ち寄る、1度だけ立ち寄る。最短の巡回経路を求めよ。というもの。

 巡回路は、1つのリングでないとならない。短い巡回路を巡り、それから突然別の町に出現し、そこから別の巡回路を巡る。2つ合わせると全ての町に立ち寄っている、というのは反則である。

 この制約条件、つまり1つのリングとなる制約条件を表現するやり方、(スリーザーリンクでの話しだが)独力では構成できなかった。しかし、実は非常に簡単なやり方があった。1つのリングであるなら、1つの矢印を除ければリングでなくなる。リングでないなら、立ち寄る順に番号を振っていける。0,1,2,…(1,2,3,…)と番号を振っていける。さらに、個別矢印の出発元の町に与えた番号は到着先の町に与えた番号よりも小さいように番号を振っていける。以下の図はその例を示している。

 

 

以下の図は、2つのリングで構成されているので、2番目のリングは「リング」のままである。

 

 

従って、2番目のリングに制約条件を満たすようには番号を与えられない。町c05の番号5は行った先の町c03の番号3よりも大きいので制約条件をみたしてない。というか、その前に町c03に何故番号3が与えられたのか理解できない。勿論、理解できない(or旨くいかない)理由は2番目の町群の矢印がリングを構成したままであるからである。

 具体的なMathProgコードを示しておこう。

 

minimize Dist: sum{i in C,j in C}d[i,j]*m[i,j];## total distance of the travel

 

s.t. M0{i in C}: m[i,i]=0;

s.t. M1{i in C}: sum{j in C} m[i,j] =1;

s.t. M2{j in C}: sum{i in C} m[i,j] =1;

 

var o1{C} integer, >=0, <=N-1;               ## o1[i]:   i-th city  o1[i]-th arrival

s.t. O0: o1["c01"]=0;

s.t. O1{ i in C, j in Cm}: o1[i]+1<=o1[j]+(N-1)*(1-m[i,j]);

#s.t. O1{ i in C, j in Cm}: o1[i]-o1[j]+N*m[i,j]<=N-1;

 

勿論、minimizeは目標。 m[i,j]は町iから町jへ行く/行かない(1/0)を示す。それと距離d[i,j]を掛けたものの総和(sum)を最小化するのが目標。

その前に制約M0,M1,M2を満たす必要がある。M0は出発元の町と到着先の町が同じであってはならないというもの。M1は、個別の町を出発する矢印があるというもの。M2は、個別の町へ到着する矢印があるというもの。制約M0, M1, M2で巡回路がリングになる保証が与えられた。ただし、複数個のリングとなる可能性がある。

1つのリングとなる保証を与えるのが制約O1である。その前の制約O0は町c01に番号0を与えているだけ。そして制約O1で個別の町に、0,1,2,…,N-1と番号を与えることになる。なお到着の町jin Cmとしてあるが、Cm= C-“c01” =c2,c3,…Cn-1である点に注意。最初の町C01に到着する矢印は考えないということである。

なお制約C12種載せている。今活きている方は、当方が普通記述するもの。素直な書き方だと、m[i,j]à o1[i]+1 =o1[j]でいいのだが、線形計画法では含意(à)を直接に指定できないので、o1[i]+1<=o1[j]+(N-1)*(1-m[i,j])とするもの。今コメント化してある方は、2つの教科書に書いてあった。当方の書き方を手で展開しただけ。ツールによっては、当方の書き方が少し遅くなる(or拒否される)可能性があるという程度の意味でしかない()

実は、1つのリングとする制約、別の書き方があった。GLPKの開発者(Andrew Makhorin)が書いているコード(exampleフォルダ内のtsp.mod)のやり方。それも示しておく。少し変えてあるが本質は同じはずである。

 

var y{C, C} >=0, <=N-1;                     ## arriving order: for avoiding sub-tour

var o0{C} >=0, <=N-1;                       ## arriving order: for dump

var o1{C} >=0, <=N-1;                       ##

s.t. CY{ i in C, j in C}: y[i,j]<=(N-1)*m[i,j];

s.t. CO{ i in Cm}: o1[i] = o0[i] +1;

 

s.t. C0{ i in C}: o0[i]= sum{j in C} y[j,i];

s.t. C1{ i in C}: o1[i]= sum{k in C} y[i,k];

 

上のコードのm[i,j]に加え、y[i,j]を導入。行く矢印がある(m[i,j]=1)場合、番号0..N-1y[i,j]に与える(制約CY)。そして、制約COでリングが1つとなる制約を与えている。この制約式内に現れた変数o0[],o1[]は、下の制約C0,C1で値が定義されているが、原本コードには現れてない。C0,C1COに展開した形が原本のコードである。ここでは、ダンプの都合でC0を導入し、ならついでということでC1も導入したもの。個別の町(C01は除く)に入ってくる矢印y[j,i]の1つだけ番号0..N-1が与えられている。その番号より1つ大きい番号を、出て行く矢印y[i,k]のどれか1つに与えるというもの。

このコード、上のコードよりも格段に速かった。10倍。さすが、GLPKの開発者ということか。GLPKの内部アルゴリズム全て分かっているわけで、遅くなることで有名な巡回セールスマン問題の推奨GLPK(or MathProg)コードを提示したということだろう。

 

[1つのリングとなる制約の記述法3] 巡回セールスマン問題の解法から

 巡回セールスマン問題は、町集合C={c01,c02,…Cn-1}が与えられ、町の間の距離d[i,j]が与えられる。町集合Cをノードとした完全グラフから、上手に辺集合を選んで、辺集合が1つの巡回路を構成する。巡回路の総移動距離が最小となるようにする。そういう問題である。従って、求解コードを書く際のわずらわしさはない。勿論、求解に時間が掛かるではあろう。

 一方、我々のスリザーリンクは、矩形領域内の格子点群が町集合に対応する。グラフの辺は隣接した格子点間に引かれる。グラフは完全グラフではありえず、(制約条件がかなりきつい)平面グラフである。ノード集合はC={c01,c02,…}などという単純な羅列とはならない。ノード集合は、[r,c] r in 1..Nr+1, c in 1..Nc+1などということになる。従って、求解コードを記述するのはうんざりする仕事になろう。しかし、計算量はたかが知れていよう。

 そういうことで、前節の[巡回セールスマン問題]で説明した「1つのリングとなる制約」の内簡単なものを採用しようと思った。しかし、それでやると、問題によっては3分(たまたま定めた許容時間)経っても終わらない問題があった。それで、うんざりする作業ではあるが、難しいほうのやり方を最終的には採用した。

 

まず、無向グラフでなく、有向グラフが必要。

 

## -------- must be one ring(directed graph ----------

        var hd0{1..Nr+1,0..Nc+1} binary; # -->  right    +-->V

        var hd1{1..Nr+1,0..Nc+1} >=0,<=1;        # <--   left    A   V

        var vd0{0..Nr+1,1..Nc+1} binary; # V     down    A<--+

        var vd1{0..Nr+1,1..Nc+1} >=0,<=1;        # A     up

s.t. C00{ r in 1..Nr+1, c in 0..Nc+1}: hd0[r,c]+hd1[r,c]=hl[r,c];

s.t. C01{ r in 0..Nr+1, c in 1..Nc+1}: vd0[r,c]+vd1[r,c]=vl[r,c];

s.t. NT2{ r in 1..Nr+1,c in 1..Nc}:hd0[r,c]<=vd1[r-1,c+1]+hd0[r,c+1]+vd0[r,c+1];

s.t. NT4{ r in 1..Nr+1,c in 1..Nc+1}: hd1[r,c]<=hd1[r,c-1]+vd1[r-1,c]+vd0[r,c]; 

s.t. NT3{ r in 1..Nr+1,c in 1..Nc+1}: vd0[r-1,c]<=hd1[r,c-1]+vd0[r,c]+hd0[r,c]; 

s.t. NT5{ r in 1..Nr+1,c in 1..Nc+1}: vd1[r,c]<=hd1[r,c-1]+vd1[r-1,c]+hd0[r,c]; 

 

変数hd1,vd1は筋からいえばbinaryだが、区間[0,1]の実数とするほうが少しは速くなるので、そうしている。

制約C00,C01は、無向辺があれば、当然有向辺があるというもの。制約NT2-5は、ノード[r,c]に到達する辺があれば出る辺もあるというもの。

 

次に、1つのリングとなる制約を指定するもの。

 

## ------- must be one ring(allocate serial number) Makhorin method------

        var o1 {1..Nr+1,1..Nc+1} >=0, <=M;       # serial number(arrival)

        var o2 {1..Nr+1,1..Nc+1} >=0, <=M;       # serial number(departure)

        var hy0{1..Nr+1,0..Nc+1} >=0, <=M;       # -->   right

        var hy1{1..Nr+1,0..Nc+1} >=0, <=M;       # <--   left

        var vy0{0..Nr+1,1..Nc+1} >=0, <=M;       # V     down

        var vy1{0..Nr+1,1..Nc+1} >=0, <=M;       # A     up

        var g  {1..Nr+1,1..Nc+1} >=0, <=1;       # nodes in the rings

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

                g[r,c]=hd0[r,c-1]+vd0[r-1,c]+hd1[r,c]+vd1[r,c];

s.t. C78{ r in 1..Nr+0, c in 1..Nc+0: ini[r,c]=9}:       # rotating right +-->V

                hd0[r,c]+vd0[r,c+1]+hd1[r+1,c]+vd1[r,c]=3;       #       A<--+

s.t. C79{ r in 1..Nr+0, c in 1..Nc+0: ini[r,c]=9}: o1[r,c+1]=1;

s.t. C80{ r in 1..Nr+1, c in 0..Nc+1}: hy0[r,c]<=M*hd0[r,c];

s.t. C81{ r in 1..Nr+1, c in 0..Nc+1}: hy1[r,c]<=M*hd1[r,c];

s.t. C82{ r in 0..Nr+1, c in 1..Nc+1}: vy0[r,c]<=M*vd0[r,c];

s.t. C83{ r in 0..Nr+1, c in 1..Nc+1}: vy1[r,c]<=M*vd1[r,c];

s.t. C84{ r in 1..Nr+0, c in 1..Nc+0: ini[r,c]<=8}: o2[r,c]=o1[r,c]+g[r,c];

s.t. C85{ c in 1..Nc+1}: o2[Nr+1,c]=o1[Nr+1,c]+g[Nr+1,c];

s.t. C86{ r in 1..Nr+1}: o2[r,Nc+1]=o1[r,Nc+1]+g[r,Nc+1];

s.t. C88{ r in 1..Nr+1, c in 1..Nc+1}:           # arival number

                o1[r,c]=hy0[r,c-1]+vy0[r-1,c]+hy1[r,c]+vy1[r,c];

s.t. C89{ r in 1..Nr+1, c in 1..Nc+1}:           # departure number

                o2[r,c]=hy1[r,c-1]+vy1[r-1,c]+hy0[r,c]+vy0[r,c];

 

 まず、制約C78などにあるIni[r,c]=9の意味。本来”3”となるべきものの1つを”9”と特別に扱うもの。

制約C70は、ノード[r,c]がリングの構成ノードかどうかを示す。

制約C78は、マス[r,c]の3辺を右回りの有向辺が回っていることを示す。リングをたどる方向を定めたことになる。

制約C79は、ノード[r,c+1]の到着番号を1とするもの。制約C78,C79は指定”9”を特別扱いしている。

制約C80-83は、有向辺に番号を与えるもの。

制約C84-86は、到着番号と出発番号の関係を示すもの。なお、リングを構成しないノードは2つの番号は共に0。

制約C88は、ノードの到着番号や出発番号と有向辺の番号との関係を与えるもの。

 

[求解時間など]

 ニコリの本(14冊目)から15個の問題を求解させた。普通のPC(3GHz,1GB)で走行。全て正解を得た。全て10*10問題。なお、最初の18*10問題も解けているが、ここでは話題にしない。

 

問題

p01

p02

p03

p04

p05

p06

p07

p08

p09

p10

p11

p12

p13

p14

p15

時間(sec)

0.6

1

0.8

6.5

1.1

15.6

2.2

6.6

2.2

9.4

1.5

3.8

0.9

1.9

2.3

3GHz

サイズ(MB)

3.7

4.2

4

4.5

4.2

4.5

4.2

4.6

3.9

4.7

4.2

4.4

3.9

4.6

4.5

 

sec、サイズはGLPKが表示したものをそのまま転載。前処理(変換)や後処理(表示)は含んでないようである。

問題によって、サイズが変わる。ただし、変わりようは小さい。

問題によって、求解時間が変わる。変わりようは大きい。また、コードを微調整するとドラスティックに時間が変わるものがある。例えばp08、上の表では6.6secとなっている。しかし、コードを少し変えると180sec(たまたま設定していた許容時間)を越えてしまう場合がある。60sec程度には直ぐになってしまう。

 この表の時間は少し掛かりすぎかとは思う。ということで、何度か微調整にトライしている。しかし、旨くはいかなかった。結局、この表の時点のコードがそれなりにいいということで、一旦、閉じたもの。

 

[高速化] 追加(130415)

 最初のコード(111018)は、速度的に少し不満であったので、今回修正を加えた。まず結果を示し、それから修正内容を示す。

 昔のコードを、今持っているPCで走行させてみた。以下の結果となった。

 

問題

p01

p02

p03

p04

p05

p06

p07

p08

p09

p10

p11

p12

p13

p14

p15

時間(sec)

0.2

0.3

0.3

2.2

0.4

5.5

0.8

2.4

0.8

3.4

0.5

1.4

0.3

0.7

0.8

i7 2.2GHz

サイズ(MB)

3.7

4.2

4

4.5

4.2

4.5

4.2

4.6

3.9

4.7

4.2

4.4

3.9

4.6

4.5

 

 周波数的には昔のPCの方が速い筈だが、キャッシュサイズの違いからだろうが、今のPCの方が3倍程度速く走ってくれている。しかし、1秒以上掛っているのが5個もある。これを軒並み1秒以下にするべく調整した。現状、以下の結果を得ている。

 

問題

p01

p02

p03

p04

p05

p06

p07

p08

p09

p10

p11

p12

p13

p14

p15

時間(sec)

0

0

0

0.4

0.1

1.9

0.7

0.9

0.4

0.4

0.2

0.3

0.2

0.3

0.4

i7 2.2GHz

サイズ(MB)

3.1

3.4

3.2

4.1

3.7

4

4.2

4.4

3.7

4.1

3.6

4.2

3.9

4.2

3.9

 

1つ(p6)、1秒以上だが、残りは1秒を切った。一応、成功したということにしよう。

 

コードの調整内容を示しておく。2つある。まず、(複数の)リングにさせる制約を書き換えた。

 

           var pa{1..Nr+1,1..Nc+1} binary;

s.t. RNG4{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c-1]<=pa[r,c];  

s.t. RNG5{ r in 1..Nr+1, c in 1..Nc+1}:vl[r-1,c]<=pa[r,c];  

s.t. RNG6{ r in 1..Nr+1, c in 1..Nc+1}:vl[r,c]  <=pa[r,c];  

s.t. RNG7{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c]  <=pa[r,c];  

s.t. RNG8{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c-1]+vl[r-1,c]+vl[r,c]+hl[r,c]=2*pa[r,c];  

 

ノードを通過するかどうかを示す変数pa[,]を導入。制約RNG8だけでも、意味論上はいいのだが、これは遅くなる。この点は、[整数計画法での解法(制約の記述法)]節で指摘しておいた。今回、制約RNG4,5,6,7を追加することで、遅くなるのを排除できた。速くなった可能性もある。

 

次に、高速化コードを導入した。2つの数字指定マスが斜め近傍に配置されている場合、

 

a,b,c,dを通過するリンクの個数を指定するもの。3,3であれば、4である。これは、有名な定理である。こういうものを指定してみた。

 

## -------- optimizations -------------------

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

           hl[r-1,c-1]+vl[r-1,c-1]+hl[r+1,c]+vl[r,c+1]=ini[r-1,c-1]+ini[r,c]-2*pa[r,c];

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

           hl[r-1,c]+vl[r-1,c+1]+hl[r+1,c-1]+vl[r,c-1]=ini[r-1,c]+ini[r,c-1]-2*pa[r,c];

 

s.t. OPT5{ r in 1..Nr, c in 1..Nc: ini[r,c]=3||ini[r,c]=9}: pa[r,c]+pa[r,c+1]+pa[r+1,c]+pa[r+1,c+1]=4;

s.t. OPT6{ r in 1..Nr, c in 1..Nc: ini[r,c]=2}: pa[r,c]+pa[r,c+1]+pa[r+1,c]+pa[r+1,c+1]>=3;

 

制約OPT1は、今説明した部分に相当する。制約OPT2は、OPT1の左右反転形に対応する。

制約OPT5は、3指定マスであれば、4隅のノード全てでリンクが通過することを示している。

制約OPT6は、2指定マスであれば、4隅のノードの内3つ以上でリンクが通過することを示している。

 

[スリザーリンクのコード]

 

以下のような起動法となる。

 

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

 

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

 

## -------------------- slither link ------------------------

        param Nr;                                # number of rows

        param Nc;                                # number of columns

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

        param M:=2*Nr*Nc -50;

## -------- must be rings ----------------------------

        var hl{1..Nr+1,0..Nc+1} binary;    # 0/1: not/link (---)

        var vl{0..Nr+1,1..Nc+1} binary;    # 0/1: not/link ( | )

        var pa{1..Nr+1,1..Nc+1} binary;

s.t. RNG0{ r in 1..Nr+1,c in 0..Nc+1 by Nc+1}: hl[r,c]=0;     # outer segments

s.t. RNG2{ r in 0..Nr+1 by Nr+1,c in 1..Nc+1}: vl[r,c]=0;     # outer segments

 

s.t. RNG4{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c-1]<=pa[r,c];  

s.t. RNG5{ r in 1..Nr+1, c in 1..Nc+1}:vl[r-1,c]<=pa[r,c];  

s.t. RNG6{ r in 1..Nr+1, c in 1..Nc+1}:vl[r,c]  <=pa[r,c];  

s.t. RNG7{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c]  <=pa[r,c];  

s.t. RNG8{ r in 1..Nr+1, c in 1..Nc+1}:hl[r,c-1]+vl[r-1,c]+vl[r,c]+hl[r,c]=2*pa[r,c];  

 

s.t. LNK3{r in 1..Nr,c in 1..Nc:ini[r,c]<=3}:    # 0,1,2,3: number of links

        hl[r,c]+vl[r,c]+hl[r+1,c]+vl[r,c+1]=ini[r,c];

s.t. LNK8{r in 1..Nr,c in 1..Nc:ini[r,c]=8}:     # 8: don't care

        hl[r,c]+vl[r,c]+hl[r+1,c]+vl[r,c+1]<=3;

s.t. LNK9{r in 1..Nr,c in 1..Nc:ini[r,c]=9}:     # 9: "3" with starting node

        hl[r,c]+vl[r,c]+hl[r+1,c]+vl[r,c+1]=3;

## -------- optimizations -------------------

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

        hl[r-1,c-1]+vl[r-1,c-1]+hl[r+1,c]+vl[r,c+1]=ini[r-1,c-1]+ini[r,c]-2*pa[r,c];

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

        hl[r-1,c]+vl[r-1,c+1]+hl[r+1,c-1]+vl[r,c-1]=ini[r-1,c]+ini[r,c-1]-2*pa[r,c];

 

s.t. OPT5{ r in 1..Nr, c in 1..Nc: ini[r,c]=3||ini[r,c]=9}: pa[r,c]+pa[r,c+1]+pa[r+1,c]+pa[r+1,c+1]=4;

s.t. OPT6{ r in 1..Nr, c in 1..Nc: ini[r,c]=2}: pa[r,c]+pa[r,c+1]+pa[r+1,c]+pa[r+1,c+1]>=3;

## -------- inhibit small ring, directly -------------------

s.t. SR0{ r in 1..Nr, c in 1..Nc-1}:     #two cells which are holzontary adjacent

        vl[r,c]+vl[r,c+2]+hl[r,c]+hl[r,c+1]+hl[r+1,c]+hl[r+1,c+1]<=5;

s.t. SR1{ r in 1..Nr-1, c in 1..Nc}:     #two cells which are vartically adjacent

        hl[r,c]+hl[r+2,c]+vl[r,c]+vl[r+1,c]+vl[r,c+1]+vl[r+1,c+1]<=5;

## ---------------- color --------------------------------

        var col{1..Nr,0..Nc} >=0,<=1;    # 0/1: white/black

s.t. COL0{ r in 1..Nr}: col[r,1]=vl[r,1];       

s.t. COL1{ r in 1..Nr, c in  2..Nc}: col[r,c]-col[r,c-1]<=vl[r,c]; # vl=0-->same

s.t. COL2{ r in 1..Nr, c in  2..Nc}: col[r,c-1]-col[r,c]<=vl[r,c]; # vl=0-->same

s.t. COL3{ r in 1..Nr, c in  2..Nc}: vl[r,c]<=col[r,c-1]+col[r,c]; # vl=1-->not

## -------- must be one ring(directed graph ----------

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

        var vd0{0..Nr+1,1..Nc+1} binary; # V     down   

        var hd1{1..Nr+1,0..Nc+1} >=0,<=1;        # <--   left   

        var vd1{0..Nr+1,1..Nc+1} >=0,<=1;        # A     up

s.t. C00{ r in 1..Nr+1, c in 0..Nc+1}: hd0[r,c]+hd1[r,c]=hl[r,c];

s.t. C01{ r in 0..Nr+1, c in 1..Nc+1}: vd0[r,c]+vd1[r,c]=vl[r,c];

s.t. NT2{ r in 1..Nr+1,c in 1..Nc}:hd0[r,c]<=vd1[r-1,c+1]+hd0[r,c+1]+vd0[r,c+1];

s.t. NT4{ r in 1..Nr+1,c in 1..Nc+1}: hd1[r,c]<=hd1[r,c-1]+vd1[r-1,c]+vd0[r,c]; 

s.t. NT3{ r in 1..Nr+1,c in 1..Nc+1}: vd0[r-1,c]<=hd1[r,c-1]+vd0[r,c]+hd0[r,c]; 

s.t. NT5{ r in 1..Nr+1,c in 1..Nc+1}: vd1[r,c]<=hd1[r,c-1]+vd1[r-1,c]+hd0[r,c]; 

## ------- must be one ring(allocate serial number) Makhorin method------

        var o1 {1..Nr+1,1..Nc+1} >=0, <=M;       # serial number(arrival)

        var o2 {1..Nr+1,1..Nc+1} >=0, <=M;       # serial number(departure)

        var hy0{1..Nr+1,0..Nc+1} >=0, <=M;       # -->   right

        var hy1{1..Nr+1,0..Nc+1} >=0, <=M;       # <--   left

        var vy0{0..Nr+1,1..Nc+1} >=0, <=M;       # V     down

        var vy1{0..Nr+1,1..Nc+1} >=0, <=M;       # A     up

        var g  {1..Nr+1,1..Nc+1} >=0, <=1;       # nodes in the rings

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

                g[r,c]=hd0[r,c-1]+vd0[r-1,c]+hd1[r,c]+vd1[r,c];

s.t. C78{ r in 1..Nr+0, c in 1..Nc+0: ini[r,c]=9}:       # rotating right +-->V

                hd0[r,c]+vd0[r,c+1]+hd1[r+1,c]+vd1[r,c]=3;       #       A<--+

s.t. C79{ r in 1..Nr+0, c in 1..Nc+0: ini[r,c]=9}: o1[r,c+1]=1;

s.t. C80{ r in 1..Nr+1, c in 0..Nc+1}: hy0[r,c]<=M*hd0[r,c];

s.t. C81{ r in 1..Nr+1, c in 0..Nc+1}: hy1[r,c]<=M*hd1[r,c];

s.t. C82{ r in 0..Nr+1, c in 1..Nc+1}: vy0[r,c]<=M*vd0[r,c];

s.t. C83{ r in 0..Nr+1, c in 1..Nc+1}: vy1[r,c]<=M*vd1[r,c];

s.t. C84{ r in 1..Nr+0, c in 1..Nc+0: ini[r,c]<=8}: o2[r,c]=o1[r,c]+g[r,c];

s.t. C85{ c in 1..Nc+1}: o2[Nr+1,c]=o1[Nr+1,c]+g[Nr+1,c];

s.t. C86{ r in 1..Nr+1}: o2[r,Nc+1]=o1[r,Nc+1]+g[r,Nc+1];

s.t. C88{ r in 1..Nr+1, c in 1..Nc+1}:      # arival number

                o1[r,c]=hy0[r,c-1]+vy0[r-1,c]+hy1[r,c]+vy1[r,c];

s.t. C89{ r in 1..Nr+1, c in 1..Nc+1}:      # departure number

                o2[r,c]=hy1[r,c-1]+vy1[r-1,c]+hy0[r,c]+vy0[r,c];

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

        var area;       # number of cells with col=1

        var length;     # nymber of links

s.t. AREA:   area  = sum{r in 1..Nr,c in 1..Nc}col[r,c];

s.t. LENGTH: length= sum{r in 1..Nr+1,c in 1..Nc+1}(hl[r,c]+vl[r,c]);

#                                        # Makhorin method      

#maximize OBJECTIVE: length-0.01*area;   # 8: 5.9 sec (6: 2.0 sec)      

maximize OBJECTIVE: length-0.1*area;     # 8: 2.3 sec (6: 3.6 sec)      

#minimize OBJECTIVE: length;     # 8: 2.3 sec (6: 3.6 sec)      

solve;

printf "length=%d\n", length;

for{r in 1..Nr}{

        for{c in 0..Nc} {

                printf "%s", if hl[r,c] then "---+" else "   +";

        } printf("\n");

                printf "%s", if vl[r,1] then "   |" else "    ";

        for{c in 1..Nc} {

                printf " %s ",if ini[r,c]=8 then "."else ini[r,c];

                printf "%s", if vl[r,c+1] then "|" else " ";

        } printf("\n");

}

        for{c in 0..Nc} {

                printf "%s", if hl[Nr+1,c] then"---+"else"   +";

        } printf("\n");

printf "--o2--\n";

for{r in 1..Nr+1}{

        for{c in 1..Nc+1} {

                printf "%4d", o2[r,c];

        } printf("\n");

}

end;

 

[問題定義ファイル]

 

ニコリのスリザーリンク第14冊内の最初の2つの問題を載せておく。"8"は" "の代用。"3"の内の1つは"9"とする。

 

#--------------- p01.txt -----------------------

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

param Nr:=10;                            # number of rows

param Nc:=10;                            # number of columns

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

        1       8 9 8 8 8 0 8 3 8 8

        2       8 0 8 1 8 3 8 0 8 2

        3       8 3 8 2 8 1 8 1 8 2

        4       8 2 8 0 8 8 8 1 8 0

        5       3 8 8 2 8 8 2 8 8 2

        6       3 8 8 0 8 8 0 8 8 3

        7       3 8 3 8 8 8 1 8 3 8

        8       2 8 3 8 2 8 0 8 3 8

        9       2 8 1 8 0 8 1 8 1 8

        10      8 8 1 8 2 8 8 8 1 8 ; /* 14-01 */

end;

#--------------- p02.txt -----------------------

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

        param Nr:=10;                            # number of rows

        param Nc:=10;                            # number of columns

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

        1       8 9 3 8 3 2 8 1 1 8

        2       3 8 8 0 8 8 1 8 8 3

        3       1 8 8 2 8 8 1 8 8 2

        4       8 3 3 8 3 3 8 0 3 8

        5       2 8 8 1 8 8 2 8 8 1

        6       2 8 8 1 8 8 2 8 8 2

        7       8 3 0 8 1 2 8 2 2 8

        8       1 8 8 1 8 8 3 8 8 3

        9       2 8 8 2 8 8 3 8 8 2

        10      8 2 3 8 3 3 8 0 2 8 ; /* 14-02 */

end;

 

 

トップページへ

 

Ads by TOK2