「充足可能性問題(3-SAT)を解く乱択アルゴリズム」 by Julia #Julia

前置き

元ネタは、結城浩氏著の「数学ガール乱択アルゴリズム」。
新しい言語を覚えるとき、慣れるために「充足可能性問題（3-SAT）を解く乱択アルゴリズム」（p.353）を実装するという癖をつけていま¹す。

ということで。前回の Egison版に引き続き。勉強開始約1ヶ月の Julia ( http://julialang.org/ ) で実装してみました²。

開発環境・動作確認環境

Mac OSX 10.9.5
Julia 0.3.5

コード

Rw3sat.jl

# Rw3sat.jl

sample(a::Array) = a[rand(1:end)]

immutable Literal
    index::Int
    not::Bool
end
literal(index, not) = Literal(index, not)

# issatisfied(l::Literal, x::BitArray{1}) = l.not ? !x[l.index] : x[l.index]
issatisfied(l::Literal, x::BitArray{1}) = l.not $ x[l.index]

immutable Clause
    literals::Array{Literal, 1}
end
clause(ls::Literal...) = Clause([ls...])

issatisfied(c::Clause, x::BitArray{1}) = any(l -> issatisfied(l, x), c.literals)

immutable CNF
    clauses::Array{Clause, 1}
end
cnf(cs::Clause...) = CNF([cs...])

issatisfied(f::CNF, x::BitArray{1}) = all(c -> issatisfied(c, x), f.clauses)

function rw3sat(f::CNF, n::Int, rr::Int)
    for r = 1:rr
        x = randbool(n) # W4
        for k = 1:3n
            if issatisfied(f, x)    # W7
                return "充足可能である"
            end
            c = sample(filter(c -> !issatisfied(c, x), f.clauses))  # W10
            idx = sample(c.literals).index  # W11
            x[idx] = !x[idx]    # W12
        end
    end
    return "おそらく充足不可能である"
end

p1 = clause(literal(1, false), literal(2, false), literal(3, false))
p2 = clause(literal(4, false), literal(2, false), literal(3, true ))
p3 = clause(literal(1, true ), literal(4, false), literal(3, false))
p4 = clause(literal(1, true ), literal(4, true ), literal(2, false))
p5 = clause(literal(4, true ), literal(2, true ), literal(3, false))
p6 = clause(literal(1, true ), literal(2, true ), literal(3, true ))
p7 = clause(literal(1, false), literal(4, true ), literal(3, true ))
p8 = clause(literal(1, false), literal(4, false), literal(2, true ))

f = cnf(p1, p2, p3, p4, p5, p6, p7, p8)
# f = cnf(p1, p2, p3, p4, p5, p6, p8)

println(rw3sat(f, 4, 3))

少し詳しく解説してみます（ツッコミ大歓迎³）。

3行目。Array から要素を乱択する関数 sample の定義。
- Julia は 1-origin。つまり配列のインデックスは 1 始まりです。
- end は、ここでは「配列の最後のインデックス」の意味になり、値としては length(a) と同じになります。
- 1:n は、Range オブジェクトを生成する書式。Ruby の 1..n、Python の range(1, n+1) 相当。
- rand(r::Range) は、Range オブジェクトが指す範囲の数値を乱択する関数。
- つまり。まとめると、a[rand(1:end)] は「配列 a のインデックス（1〜length(a)）を1つ乱択して、その位置の要素を返す」という意味になるわけです。
5〜8行目。変更不可能（immutable）な型 Literal の定義。index（Int型、変数のインデックス）と not（Bool型、true ならそのまま、false なら真偽値反転）の2つのメンバからなる型として定義しています。
- ここでは immutable キーワードを利用していますが、変更可能（mutable）な型を定義する場合は type キーワードを使用します。
9行目。Literal 型の外部コンストラクタ定義。ここでは受け取った引数をデフォルトコンストラクタに引き渡しているだけです。
11〜12行目。issatisfied の定義。第1引数 l（Literal型）と x（BitArray{1}型）を受け取って、その Literal が満たされているかどうかを返します。
- 11行目が意味的に妥当な丁寧な実装。12行目はその shorthand。$ は Julia の排他的論理和演算子。つまり「l.not と x[l.index] の真偽値が一致しないときのみ true」というワケ。
- BitArray は、1要素1ビットからなる配列。要素は Bool 値、だけど Array{Bool} とは別物（←をパックしたもの）。
14〜16行目。immutable 型 Clause の定義。Literal の配列をメンバとして持つ。
17行目。Clause の外部コンストラクタ定義。
- 左辺：仮引数の ls::Literal... は、所謂「可変長引数」の書式。ls は、全ての要素が Literal型となる tuple になります。Python や Ruby の引数 *ls という書式と同様（Python では同様に tuple として受け取れる）。
- 右辺：デフォルトコンストラクタに渡している実引数の [ls...] という書式は、iterable な変数 ls を展開して配列のコンストラクタに渡しているカタチ。Rubyでも [*ls] と書けますね。
- まとめると、c = clause(l1, l2, l3) と書くと、c = Clause([l1, l2, l3]) として処理され、c.literals == [l1, l2, l3] となる、ということです。
19行目。issatisfied の定義（多重定義）。第1引数が Clause で、第2引数は先ほどと同じ BitArray{1}。
メンバ literals のいずれか1つでも issatisfied ならば true となります。
- この「関数多重定義」が Julia の（たぶん最大の）特長。
  Julia では、「メソッドは関数に属する」という思想。「新しい型を定義して、その型を引数に取る関数を（多重）定義する（＝メソッドを定義する）」、というプログラミングスタイル（なんだと思っている）。
  例えば今回定義していませんが、新しく定義した型の文字列表現定義（print 関数等に渡したときに表示される形式の定義）は、(Base.)show という関数を多重定義します。この行為が、Rubyでいう to_s メソッドの定義、Python なら __str__() の定義に相当します。
- l -> 〜 は lambda 式（無名関数）の定義。Python なら lambda l: 〜、Ruby なら ->l {〜} に相当。
21〜23行目。immutable 型 CNF の定義。Clause の配列をメンバとして持つ。
24行目。CNF の外部コンストラクタ定義。Clause の時と同様なので詳細略。
- f = cnf(c1, c2, c3) と書くと、f = CNF([c1, c2, c3]) として処理され、f.clauses == [c1, c2, c3] となります。
26行目。issatisfied の定義（多重定義）。第1引数が CNF で、第2引数は先ほどと同じ BitArray{1}。
メンバ clauses のすべてが issatisfied ならば true となります。
28〜41行目。3-SAT の本体。Ruby や Python のようなループ制御・条件分岐の構文が存在するので、やってること自体はほぼ書籍の仮想コードそのままです。
今回も対応する仮想コードの行番号を何ヶ所か記述しておいたので、書籍をお持ちの方はそれと対比して見ていただけると良いと思います。
それ以外の、Julia 固有の部分について詳説↓：
- 29行目・31行目。for r = 1:rr という書き方は（どうやら）MATLAB/Octave 由来。R や Python ライクに for r in 1:rr と書くこともできます（同じ意味になります）。
- 30行目。randbool(n::Integer) は、要素数 n の BitArray を返す関数。各要素は false か true かが乱択される仕組み。
- 31行目。3n は 3 * n の shorthand。数字と識別子が間に空白を挟まずに続けて記述されていたら積と解釈される仕様。これを利用して例えば方程式が 3x^2 - 2x + 1 のように数学的に直感的に書けます（^ はべき乗演算子）。
43〜50行目。各 Clause の定義。
52(53)行目。CNF の定義。
55行目。SAT を実行して結果を表示。
- 52行目有効（53行目コメントアウト）なら「おそらく充足不可能である」、53行目が有効なら「充足可能である」と表示されます。

Julia の技術情報的な内容以外（3-SAT そのものについて等）は、書籍や、私や他の方の過去の実装例（↓参照）を参考にしてください。

参考

私が過去に実装した他言語での実装例

Ruby版: https://gist.github.com/960956
JavaScript版: https://gist.github.com/966408
CoffeScript版: https://gist.github.com/1000913
R版: https://gist.github.com/2811277
Python版: https://gist.github.com/2818831
Haxe版: https://gist.github.com/4148331
Egison版: http://qiita.com/antimon2/items/c55a31a1bbd1f98360f0

その他の実装例・関連記事

Julia の技術情報

Julia Documentation（公式ドキュメントけっこう充実）
Julia Advent Calendar 2014 の各記事（これがなかったら1ヶ月でここまで辿り着けなかったと思います多謝）（一括ご容赦m(_ _)m）
その他 julialang.org 内参照。

以上、前回記事からのコピペ。 ↩
実際に実装したのは少し前、この記事を書き始めたのは昨日（2015/02/01）なのですが、事情があって寝落ちしました…。 ↩
まだまともに触りだして1ヶ月なのでいろいろ理解が間違っているところとか解釈がおかしいところとかはあると思います。指摘があり次第順次この記事を修正していく予定です。もちろん「ここはもっとこう書いた方が良いよ」というアドバイスも大歓迎です。 ↩