Wu らによる差分検出の O(NP) アルゴリズム実装シリーズ ③Goによる実装

はじめに

これは、3回にわけて Wu らによる O(NP) アルゴリズムを解説、実装するシリーズです。

• ①仕組み、考え方
• ②探索方法
• ③Goによる実装 <- 今回

の3部に分かれています。

今回は、O(NP) アルゴリズムを実際に実装していきます。
前回(①仕組み、考え方)、前々回(②探索方法)でアルゴリズムの考え方、探索方法はわかっているので、今回は実装上の注意点を説明しつつ、実装していきます。

元論文は An O(NP) Sequence Comparison Algorithm by described by Sun Wu, Udi Manber and Gene Myers です。興味のある方はそちらも目を通すと面白いと思います。

最終的な実装を確認したい方は https://github.com/convto/onp に置いてあるのでそちらをごらんください。

これは Makuake Development Team Advent Calendar 2019 - Adventar の21日目の投稿です(遅くなりました)

前回のおさらい

②探索方法 で最終的に整理した擬似コードを振り返ります。
今回は、この擬似コードをベースに実装を始めます

dummy.code

a = 入力A
b = 入力B

for (p = 0; p <= M, p++) {
    k ∈ [-p, Δ+1] {
        fp[k][p] = (k, max(fp(k-1, p) + 1, fp(k+1, p-1)))
    }
    k = Δ {
        fp[k][p] = snake(k, max(k-1, p) + 1, max(k+1, p))
    }
    k ∈ [Δ+1, Δ+p] {
        snake(k, max(fp(k-1, p-1) + 1, fp(k+1, p)))
    }
    if (fp[Δ][p] == N) {
        return Δ + 2*P // Dの値を返す
    }
}

snake(k, y) int {
    x = y - k
	while (a[x] == b[y]) {
		x++
		y++
	}
	return y
}

Goによる実装

さきほどの擬似コードのような処理を実装すれば最短編集距離が導けます。
この手順ではSESやLCSなどは求めていないですが、それぞれ

• fp(k-1, p) と fp(k+1, p) のどちらを採用したのか覚えておいて 点(M,N) にたどり着いたときに順番通り出力すれば SES が求められる
• snake で斜めに進めた(a[x] と b[y] が一致した)かどうかを覚えておいて 点(M,N) にたどり着いたときに一致した文字列をつなげれば LCS が求められる

などの方法でサクッと求められるので、今回の説明ではこれらを求める部分については省略します。

では実装です。早速、先程の擬似コードのかんたんに実装できる部分についてコード化していきます。
(一部擬似コードを残しているので実行できません)

package onp

func max(x, y int) int {
	if x > y {
		return x
	}
	return y
}

func onp(a, b []rune) int {
    m := len(a)
    n := len(b)
    // ①
	for p := 0; ; p++ {
		k ∈ [-p, Δ+1] {
    	    fp[k][p] = (k, max(a, b, m, n, fp(k-1, p) + 1, fp(k+1, p-1)))
    	}
    	k = Δ {
    	    fp[k][p] = snake(a, b, m, n, k, max(k-1, p) + 1, max(k+1, p))
    	}
    	k ∈ [Δ+1, Δ+p] {
    	    snake(k, max(a, b, m, n, fp(k-1, p-1) + 1, fp(k+1, p)))
    	}
    	if (fp[Δ][p] && kp[k][p] == N) {
    	    return Δ + 2*P // Dの値を返す
    	}
    }
}

// k, y 以外にも処理のために a,b,m,n が必要
func snake(a, b []rune, m, n, k, y int) int {
	x := y - k
    // ②
	for x < m && y < n && a[x] == b[y] {
		x++
		y++
	}
	return y
}

注意としては
① P は M 以下ですが、Mまでに 点(M,N) に到着するのであえて制御する必要はないです
② 入力A,Bの長さM,Nを越えた値をx,yが取るとエディットグラフをはみ出てしまうので、 x < m または y < n のときは右下に進みません。

k の探索のコード化

後は k についての探索をコード化すれば良さそうです。
まず、snake の結果を入れる k 上の最遠点 fp(k,p) を表現する表現する配列として fp []int を置きます。
ここで、 fp(k,p) から値を取り出すには k, p が必要なので、 fp [][]int などの二次元配列でないと表現できないのでは？と思った方もいらっしゃると思います。

これは、実装ですこし工夫することで　fp(k,p) の値を k だけで取り出すことができるため一次元配列としています。

var fp []int // 初期化は別途整理
for p := 0; ; p++ {
	// ①
	for k := -p; k <= delta-1; k++ {
		fp[k] = snake(k, max(fp[k-1]+1, fp[k+1]))
	}
	// ②
	for k := delta + p; k >= delta+1; k-- {
		fp[k] = snake(k, max(fp[k-1]+1, fp[k+1]))
	}
	// ③
	fp[delta] = snake(delta, max(fp[delta-1]+1, fp[delta+1]))
	if fp[delta] == n {
		return delta + 2*p
	}
}

ある値 p のときの k の配列、というのはこのように1次元配列で表現できます。
このとき, k ∈ [-p, Δ-1] のときと k ∈ [Δ+1, Δ+p] のときとでループのの方向が違うことに注目してください。
なぜループの方向を変えているかというと、

          snake(k, max(fp(k-1, p) + 1, fp(k+1, p-1)))  (k ∈ [-p, Δ-1])
fp(k,p) { snake(k, max(k-1, p) + 1, max(k+1, p))       (k = Δ)
          snake(k, max(fp(k-1, p-1) + 1, fp(k+1, p)))  (k ∈ [Δ+1, Δ+p])

より

• k ∈ [-p, Δ-1] では k+1 方向で一つ前のpのときの最遠点が知りたい
• k ∈ [Δ+1, Δ+p] では k-1 方向で一つ前のpのときの最遠点が知りたい

からです。

一つ前のpのときの最遠点をループの方向で表現することで、fp(k,p) を fp []int の一次元配列で表現することができます。

①では snake(k, max(fp(k-1, p) + 1, fp(k+1, p-1))) (k ∈ [-p, Δ-1]) より、 k+1 の最遠点 fp[k+1] は一つ前のpのループでの値を参照したいです。
k をインクリメントしてループさせることにより、fp[k+1] はまだ求められていないので一つ前のpループでの値を参照できます。

同じように②では snake(k, max(fp(k-1, p-1) + 1, fp(k+1, p))) (k ∈ [Δ+1, Δ+p]) より k-1 の最遠点 fp[k-1] は一つ前のpのループでの値を参照したいです。
k をデクリメントしてループさせることにより、 fp[k-1] はまだ求められていないので一つ前のpループでの値を参照できます。

次に③のときですが、
これは snake(k, max(k-1, p) + 1, max(k+1, p)) (k = Δ) より p-1 を参照しないのでループの方向などを考慮する必要はないです。(そもそもループしませんが)
一部①、②の結果を参照するのでこれらより後に実行する必要があります。

添字のオフセット設定、初期化

このままだとfpを取り出すときの添字がマイナスになるケースがあります。
Goの場合はスライスにマイナスの添字でアクセスするとpanicが起きるので、うまくマイナス分を相殺するためにオフセットを取ります。

fp []int がアクセスされうる最も大きいマイナス値は -(M+1) です。
これは入力AとBの文字列が完全に一致しないためBの文字長文すべて削除するなどのケース、P=M となる場合です。

なぜ +1 があるかというと、先程のコードの①のときに fp[k-1] で k=-P=-M のときにfpの添字が fp[(-M)-1] となるためです。
それらを打ち消すためには M+1 のオフセットがあればマイナス添字へのアクセスを避けられるので、オフセット値を M+1 とします

offset := M + 1

次にfpのサイズを初期化します。fpの添字はkなので、kの取りうる最大の値と最小の値の範囲をすべて格納できるようにサイズ確保します。

k が取りうる一番小さい値は先程説明したように -(M+1) です。

k が取りうる一番大きい値は、入力Bが空文字で空文字で全文字挿入するなどのケース、k=P=N となる場合です。
先程の実装の②のとき fp[k+1] を参照するので fp[(N)+1]となる場合が添字の値が一番大きくなります。

一番小さい値はマイナス値のため、先程説明したようにオフセットが加算されます。
つまり取りうる添字で一番大きい値は

offset + (N+1) = (M+1) + (N+1)
               = M+N+2

となります。

添字は0から始まるので、この取りうる最大の添字の値に添字0のときの分で1足してあげればサイズは求められます。
つまりfpの初期化は以下のように書けます

fp := make([]int, m+n+3)

これまでの実装をまとめると以下のようになります

package onp

func max(x, y int) int {
	if x > y {
		return x
	}
	return y
}

// len(b) >= len(a)
func onp(a, b []rune) int {
	m, n := len(a), len(b)
	offset := m + 1
	delta := n - m
	fp := make([]int, m+n+3)

	for p := 0; ; p++ {
		for k := -p; k <= delta-1; k++ {
			fp[k+offset] = snake(a, b, m, n, k, max(fp[k-1+offset]+1, fp[k+1+offset]))
		}
		for k := delta + p; k >= delta+1; k-- {
			fp[k+offset] = snake(a, b, m, n, k, max(fp[k-1+offset]+1, fp[k+1+offset]))
		}
		fp[delta+offset] = snake(a, b, m, n, delta, max(fp[delta-1+offset]+1, fp[delta+1+offset]))
		if fp[delta+offset] == n {
			return delta + 2*p
		}
	}
}

// k, y 以外にも処理のために a,b,m,n が必要
func snake(a, b []rune, m, n, k, y int) int {
	x := y - k
	for x < m && y < n && a[x] == b[y] {
		x++
		y++
	}
	return y
}

ここまでくると一旦実行はできるので、パッケージをmainにしてmain関数を追加すれば編集距離が求められます
https://play.golang.org/p/xQj0TDEHsn1

値の初期化

さて、ここで snake にわたすy値を取得するために max で k-1 と k+1 の最遠点を比較する部分ですが、
まだ探索していない値には、Goの場合はゼロ値である0が入っています。
探索した結果、あるkの最遠点がのy値が0となるパターンは発生しうります。

つまり

• 実際に最遠点が0なのか
• 初期値の0なのか

がわからない場合があります。

編集距離を求めるだけならyの値さえわかればいいので、値が同値なら k-1 と k+1 のどちらの最遠点を採用したかはわからなくても問題になりません。
ですが、SESなどを取得しようとすると、どちらから採用したのかわからないと挿入だったのか削除だったのかわからなくなってしまいます。

そのため、論文では p<0 または k ∉ [-p, Δ+p] のとき fp(k,p) に -1 を代入する方法を紹介しています。
-1 と 0 なら 0 のほうが大きいので、最遠点が0になるパターンでも問題なく挿入と削除を区別できます。

p の値は探索しないとわからない値で、pは0から始まるので fp []int のすべての要素を -1 で初期化すればよいです。

fp := make([]int, d.m+d.n+3)
for i := range fp {
	fp[i] = -1
}

これを先程のコードに当てはめつつ、snakeの引数が少し多い部分の整理や M <= N の保証、関数の公開など対応すると以下のようになります。

package onp

type Diff struct {
	a, b []rune
	m, n int // m <= n, m is a length, n is b length
}

func max(x, y int) int {
	if x > y {
		return x
	}
	return y
}

// NewDiff creates Diff so that m <= n.
func NewDiff(a, b []rune) Diff {
	m, n := len(a), len(b)
	if m > n {
		return Diff{a: b, b: a, m: n, n: m}
	}
	return Diff{a: a, b: b, m: m, n: n}
}

func (d Diff) ONP() int {
	offset := d.m + 1
	delta := d.n - d.m
	fp := make([]int, d.m+d.n+3)
	for i := range fp {
		fp[i] = -1
	}
	for p := 0; ; p++ {
		// -p <= k <= delta - 1
		for k := -p; k <= delta-1; k++ {
			fp[k+offset] = d.snake(k, max(fp[k-1+offset]+1, fp[k+1+offset]))
		}
		// d.delta + 1 <= k <= d.delta + p
		for k := delta + p; k >= delta+1; k-- {
			fp[k+offset] = d.snake(k, max(fp[k-1+offset]+1, fp[k+1+offset]))
		}
		// d.delta == k
		fp[delta+offset] = d.snake(delta, max(fp[delta-1+offset]+1, fp[delta+1+offset]))
		if fp[delta+offset] == d.n {
			return delta + 2*p
		}
	}
}

// snake denote the y-coordinate of the furthest point on diagonal k
func (d Diff) snake(k, y int) int {
	x := y - k
	for x < d.m && y < d.n && d.a[x] == d.b[y] {
		x++
		y++
	}
	return y
}

この最終的なコードは https://github.com/convto/onp においてあります。

これで実装できました。
ループの方向を工夫したり、要素列のサイズを計算したり、すこし複雑な操作もありましたが比較的素直にアルゴリズムどうりに実装できたと思います。

さいごに

Wu らによる O(NP) アルゴリズムを紹介しました。
全探索に比べて比べて計算量が少ない理由や、探索方法についても概観がつかめたと思います。

筆者自身はまともに論文を読んだのはこれが初めてで、最初は主張している内容も実装もどちらも全く理解できなかったです。
しっかり理解できるまで一生懸命論文や実装を見て、手を動かしてみて...など繰り返してようやくわかってきたところです。

主張やそれが成立する理由がわかったときにはすごくスッキリしました。理解できたのがすごく嬉しかったことを覚えています。
今回は一人でも多くの人にそのスッキリ感を、できる限りわかりやすく伝えよう！と思い筆をとった次第です。
どこまでわかりやすく書けたかは正直わかりませんが、この文章がいつか誰かの助けになれば幸いです。

また、この内容は同時進行で社内の勉強会勉強会でも解説していて、そのような学びの機会が社内にあるのもとてもありがたいです！

参考文献

• https://www.sciencedirect.com/science/article/pii/002001909090035V 元論文
• https://github.com/cubicdaiya/gonp 今回は実装しませんでしたが、LCSやSESを求めるときの処理などで参考にしました。わかりやすかったです...！