Edited at

最適化における Julia - JuMP で最長しりとり?

More than 1 year has passed since last update.

 Julia の最適化パッケージ JuMP を試してみる連載。

 Qiita記事 「最長しりとりを組合せ最適で解く」 @Tsutomu-KKE@github さんによる Python - PuLPプログラムを、JuMP に移植してみます。


最適化プログラム


using JuMP
using Clp

function siritori(kw)
n=length(kw)
m=Model()
@variable(m, x[1:n,1:n], Bin)
@variable(m, y[1:n], Bin)
@objective(m, Max, sum(x[i,j] for i=1:n, j=1:n)) # (0) なるべく繋げる
@constraint(m, [i=1:n], x[i,i] == 0) # (A)
for i=1:n, j=1:n
if kw[i][end] != kw[j][1]
@constraint(m, x[i,j] == 0) # (1)
end
end
@expression(m, cou[i=1:n], sum(x[i,j] for j=1:n) )
@expression(m, cin[i=1:n], sum(x[j,i] for j=1:n) )
@constraint(m, [i=1:n], cou[i] <= 1) # (3) kw_i から出る数は1以下
@constraint(m, [i=1:n], cin[i] <= 1) # (4) kw_i へ入る数は1以下
@constraint(m, [i=1:n], cou[i] <= cin[i] + y[i] ) # (5) yに関する制約
@constraint(m, sum(y[i] for i=1:n) == 1) # (6) 先頭は1つだけ
@time status=solve(m)
println(status)
ys=getvalue(y)
xs=getvalue(x)
return m, status, ys, xs
end

 関数 siritori(kw) は、文字列(単語)の配列 kw に対して、最適化を用いて、最長しりとりを求めようとするものです。

 定式化は参照記事 をご覧ください。プログラム中の (0)-(6) は、参照記事の条件に対応します。

 実装上の大きな違いは、定義する変数です。参照記事では「しりとり」が可能な単語間の向き付きリンク ($i \rightarrow j$ ) に対してのみ変数 $x_{i,j}$ を生成するのに対して、ここでは、変数 $x_{i,j}$ を全て生成した上で、「しりとり」が不可能なリンクに対して制約条件で $x_{i,j} = 0$ を課しています。

 更に、自分自身へのリンク ($i \rightarrow i$) を、明示的に禁止しました。$x_{i,i} = 0$。プログラムコメント: 制約(A)


簡単な例

簡単な例で実験します。

julia> kw1 = [ "cde", "txy", "abc", "rst", "rsr", "pqr"  ]

6-element Array{String,1}:
"cde"
"txy"
"abc"
"rst"
"rsr"
"pqr"

julia> mm1, stat1, yo1, xo1 = siritori(kw1)
0.015661 seconds (99 allocations: 24.438 KB)
Optimal
(Maximization problem with:
* 55 linear constraints
* 42 variables: 42 binary
Solver is default solver,:Optimal,[0.0,0.0,0.0,0.0,0.0,1.0],
[0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; ; 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0])

julia> getobjectivevalue(mm1)
3.0

julia> yo1
6-element Array{Float64,1}:
0.0
0.0
0.0
0.0
0.0
1.0

julia> xo1
6×6 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0 1.0 0.0

 検出された「しりとり」(リンク)を印字する関数を書きましょう。

function print_links(kw,xs)

count=0
n=size(xs,1)
for is=1:n
ie=findfirst(xs[is,:])
if (ie > 0)
println( "$(is) $(kw[is]) -> $(ie) $(kw[ie])" )
count += 1
end
end
println("count=$(count)")
end

 $x_{i,j} = 1$ なら、$i \rightarrow j$ の「しりとり」(リンク)が成立しています。各行 $i$ のスライス x[i,:] から 1 が入っている列番号 $j$ を探索します。探索には、関数 findfirst(v) を用います。これは、0 以外の数が最初に見つかった配列添字を返します。見つからなければ戻り値は 0 です。

 (式1) || (式2) は、



if (not (式1))

(式2)

end

と等価です。 (式1) || (式2) の論理式の真偽は、(式1) が偽なら、(式2) を調べないと分からないからです。 short circuit evaluation といいます。

 文字列中に "$(式)" と書くと、式の評価結果で置き換えられます。(Ruby の "#{式}" と同じですね。)

 実行してみましょう。

julia> print_links(kw1,xo1)

4 rst -> 2 txy
5 rsr -> 4 rst
6 pqr -> 5 rsr
count=3

 「最長しりとり」を印字する関数を書きます。「しりとり」の最初の単語の番号は、配列 y が非零となる添字です。

function print_siri(kw,ys,xs)

is=findfirst(ys)
is > 0 || return
print( "$(is) $(kw[is])")
count=0
while true
is=findfirst(xs[is,:])
is > 0 || break
print( " ->\n$(is) $(kw[is])" )
count += 1
end
println("\ncount=$(count)")
end

 では実行。「最長しりとり」が検出されていますね。

julia> print_siri(kw1,yo1,xo1)

6 pqr ->
5 rsr ->
4 rst ->
2 txy
count=3


少し長い例。

 リンク先にあるように、C++のキーワード 84個で「しりとり」をしてみます。

 まずキーワードを作ります。

julia> cppkw = """

alignas,alignof,and,and_eq,asm,auto,bitand,bitor,bool,break,case,
catch,char,char16_t,char32_t,class,compl,const,constexpr,const_cast,
continue,decltype,default,delete,do,double,dynamic_cast,else,enum,
explicit,export,extern,false,float,for,friend,goto,if,inline,int,long,
mutable,namespace,new,noexcept,not,not_eq,nullptr,operator,or,or_eq,
private,protected,public,register,reinterpret_cast,return,short,
signed,sizeof,static,static_assert,static_cast,struct,switch,template
this,thread_local,throw,true,try,typedef,typeid,typename,union,
unsigned,using,virtual,void,volatile,wchar_t,while,xor,xor_eq"""

"alignas,alignof,and,and_eq,asm,auto,bitand,bitor,bool,break,case,\ncatch,char,char16_t,char32_t,class,compl,const,constexpr,const_cast,\ncontinue,decltype,default,delete,do,double,dynamic_cast,else,enum,\nexplicit,export,extern,false,float,for,friend,goto,if,inline,int,long,\nmutable,namespace,new,noexcept,not,not_eq,nullptr,operator,or,or_eq,\nprivate,protected,public,register,reinterpret_cast,return,short,\nsigned,sizeof,static,static_assert,static_cast,struct,switch,template\nthis,thread_local,throw,true,try,typedef,typeid,typename,union,\nunsigned,using,virtual,void,volatile,wchar_t,while,xor,xor_eq"

julia> ckw=split(cppkw, r"[,\n]+")
84-element Array{SubString{String},1}:
"alignas"
"alignof"
"and"
"and_eq"
"asm"
"auto"
"bitand"
"bitor"
"bool"
"break"

"unsigned"
"using"
"virtual"
"void"
"volatile"
"wchar_t"
"while"
"xor"
"xor_eq"

 """ で始める文字列は改行できます。 """ で終了します。改行文字 \n が入ることに注意しましょう。

  split(s, ch) は、文字列 s を、文字や正規表現 ch で区切って、文字列配列を作ります。ここでは、カンマ , または改行文字で区切りました。

 実行してみましょう。

julia> cm, c, cy, cx = siritori(ckw)

0.148403 seconds (134 allocations: 3.093 MB)
Optimal
(Maximization problem with:
* 6975 linear constraints
* 7140 variables: 7140 binary
Solver is default solver,:Optimal,[1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 … 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

julia> print_links(ckw,cx)

1 alignas -> 59 signed
16 class -> 60 sizeof
23 default -> 73 typeid
25 do -> 50 or
27 dynamic_cast -> 72 typedef
28 else -> 30 explicit
29 enum -> 42 mutable
30 explicit -> 74 typename
31 export -> 67 this
32 extern -> 44 new
34 float -> 66 template
36 friend -> 27 dynamic_cast
37 goto -> 49 operator
41 long -> 37 goto
42 mutable -> 28 else
44 new -> 82 while
46 not -> 68 thread_local
49 operator -> 57 return
50 or -> 55 register
55 register -> 56 reinterpret_cast
56 reinterpret_cast -> 69 throw
57 return -> 46 not
59 signed -> 23 default
60 sizeof -> 36 friend
61 static -> 16 class
67 this -> 61 static
68 thread_local -> 41 long
69 throw -> 81 wchar_t
70 true -> 31 export
72 typedef -> 34 float
73 typeid -> 25 do
74 typename -> 32 extern
81 wchar_t -> 70 true
82 while -> 29 enum
count=34

julia> getobjectivevalue(cm)
34.0

 34個のリンクが得られました。

 最長しりとりを印字しましょう。

julia> print_siri(ckw,cy,cx)

1 alignas ->
59 signed ->
23 default ->
73 typeid ->
25 do ->
50 or ->
55 register ->
56 reinterpret_cast ->
69 throw ->
81 wchar_t ->
70 true ->
31 export ->
67 this ->
61 static ->
16 class ->
60 sizeof ->
36 friend ->
27 dynamic_cast ->
72 typedef ->
34 float ->
66 template
count=20

あれれ、20回目で途切れてしまいました。

print_links の結果から、手で「しりとり」してみます。

1 alignas -> 59 signed

59 signed -> 23 default
23 default -> 73 typeid
73 typeid -> 25 do
25 do -> 50 or
50 or -> 55 register
55 register -> 56 reinterpret_cast
56 reinterpret_cast -> 69 throw
69 throw -> 81 wchar_t
81 wchar_t -> 70 true
70 true -> 31 export
31 export -> 67 this
67 this -> 61 static
61 static -> 16 class
16 class -> 60 sizeof
27 dynamic_cast -> 72 typedef
28 else -> 30 explicit
29 enum -> 42 mutable
30 explicit -> 74 typename
32 extern -> 44 new
34 float -> 66 template

41 long -> 37 goto
37 goto -> 49 operator
49 operator -> 57 return
57 return -> 46 not
46 not -> 68 thread_local
68 thread_local -> 41 long

44 new -> 82 while
82 while -> 29 enum

36 friend -> 27 dynamic_cast

42 mutable -> 28 else

60 sizeof -> 36 friend

72 typedef -> 34 float

74 typename -> 32 extern

リンクは、以下に分類されました。



  • 1 から 66


  • 41 からの循環。


  • 448229

  • 続かない「しりとり」 5個

 振り返ってみれば、今回の定式化は、「"しりとり"の総数を最大にすること」が目的であり、「しりとりをなるべく連結する」という条件は入っていません。$x_{i,j}$ の組合せは有限ですから、全ての場合を試行すれば 「最長しりとり」を見つけられるはずですが、今回の目的関数では適当ではなさそうです。


リンク

Julia の最適化パッケージ JuMP に関する記事です。

- 「最適化における Julia : JuMP事始め」

- 「同: JuMP で N-Queen を解く」

- 「同: - JuMP で色々解く」 (今後も追記予定)