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

　Julia の最適化パッケージ JuMP を試してみる連載。
　Qiita記事 [「最長しりとりを組合せ最適で解く」](http://qiita.com/Tsutomu-KKE@github/items/ba2dedb5795cae36f8a1) @Tsutomu-KKE@github さんによる Python - PuLPプログラムを、JuMP に移植してみます。

# 最適化プログラム
```Julia

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` に対して、最適化を用いて、最長しりとりを求めようとするものです。
　定式化は[参照記事](http://qiita.com/Tsutomu-KKE@github/items/ba2dedb5795cae36f8a1) をご覧ください。プログラム中の (0)-(6) は、参照記事の条件に対応します。
　実装上の大きな違いは、定義する変数です。参照記事では「しりとり」が可能な単語間の向き付きリンク ($i \rightarrow j$ ) に対してのみ変数 $x_{i,j}$ を生成するのに対して、ここでは、変数 $x_{i,j}$ を全て生成した上で、「しりとり」が不可能なリンクに対して制約条件で $x_{i,j} = 0$ を課しています。
　更に、自分自身へのリンク ($i \rightarrow i$) を、明示的に禁止しました。$x_{i,i} = 0$。プログラムコメント: 制約(A)

# 簡単な例

 簡単な例で実験します。

```Julia
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
```

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

```Julia
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)`](http://docs.julialang.org/en/stable/stdlib/arrays/#Base.findfirst) を用います。これは、`0` 以外の数が最初に見つかった配列添字を返します。見つからなければ戻り値は `0` です。
　`(式1) || (式2)` は、
```
if (not (式1))
  (式2)
end
```

と等価です。 `(式1) || (式2)` の論理式の真偽は、`(式1)` が偽なら、`(式2)` を調べないと分からないからです。 [short circuit evaluation](http://docs.julialang.org/en/stable/manual/control-flow/#short-circuit-evaluation) といいます。
　文字列中に `"$(式)"` と書くと、式の評価結果で置き換えられます。(Ruby の `"#{式}"`　と同じですね。)
　実行してみましょう。

```Julia
julia> print_links(kw1,xo1)
4 rst -> 2 txy
5 rsr -> 4 rst
6 pqr -> 5 rsr
count=3
```

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


```Julia
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
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)`](http://docs.julialang.org/en/stable/stdlib/strings/#Base.split) は、文字列 `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 
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` からの循環。
- `44` → `82` → `29` 
- 続かない「しりとり」 5個

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

# リンク

Julia の最適化パッケージ JuMP に関する記事です。
- [「最適化における Julia : JuMP事始め」](http://qiita.com/tenfu2tea/items/499ff93897f25a0256c0) 
- [「同: JuMP で N-Queen を解く」](http://qiita.com/tenfu2tea/items/b1ad78cc9f0942fdbc93)
- [「同: - JuMP で色々解く」](http://qiita.com/tenfu2tea/items/29cb927f1e785047532c) (今後も追記予定)