枕を欹てて聴く

香炉峰の雪は簾を撥げて看る

Diff algorithm


id:smoking186 さんの指摘を受け, First Authorの名前などを付加しました.
どうもです.

記事内のcodeは最適化などを施しておらず, 冗長に, 定義どおりに書いています.
ifがまとめられたりとかしますが, そのあたりはご容赦を...

Rubyでlevenshtein距離を見て以来, 個人的にdiffブームが来ていた.
計算量O(ND) / O(NP)のalgorithmなどがあるのは知っていたが, 論文(英語)および, 解説のみ, またはソースコードのみなど分かれているものが多く, algorithmに疎い自分には理解するのに大変時間がかかってしまった. しかしやっとわかったので, 解説+JS実装してみる. 解説とソースコードがセットだと, 多少はわかりやすくなるかと...
自分は正直これくらい細かく言われないとすぐにはわかんない人なので(the O(ND)だけで2日かかった), かなり冗長に書いています.

二次元配列全走査によるLevenshtein距離の求め方及び, それを使ったdiff

Levenshtein距離

Levenshtein距離は, ある文字列Aからある文字列Bへ変換する際にどれだけの編集回数が必要かというものを数値化したものです. これが小さいと, 編集回数が近い, つまり文字列の類似性が高いということがわかります.
例えば,
string を strength に変換する際には

| s
| t
| r
- i
+ e
| n
| g
+ t
+ h

より, -が1回, +が3回で合計4が編集距離となります.
これを求めるやり方に, EditGraphというものを使います.
f:id:Constellation:20091021161357p:image
上のように, 横列に文字列A, 縦列に文字列Bをおきます. そして,

  • 縦に降りる処理 = 文字を+する
  • 横にずれる処理 = 文字を-する

という風にとると, 編集の道筋を取ることができます.

そして次に, 文字が同じ線が交差する点を斜線で結びます.
f:id:Constellation:20091021161358p:image
この斜線部分は文字が同じであり, ここを通る処理 = 文字列の改変しない という風に考えることができます.
そして,

  • 文字を+する処理はcost 1
  • 文字を-する処理はcost 1
  • そのままにする処理をcost 0

と考えると, たどった経路のcost = 編集距離とすることができます.

このEditGraphの右下へたどり着いたとき, 編集が終了するということは, 編集距離を求める行為は, いかに最小costで右下にたどり着くか. もしくは言い換えれば, いかにcost 0である斜線部分を通りまくるって右下にたどり着くかという問題にすることができます.

では実際に求めて見ます.
f:id:Constellation:20091021161356p:image
上の図のように, ある点に行く点は,

  • 上から降りる
  • 横からずれる

の2つが基本的には存在し, 文字が同じときだけボーナス経路として

  • 斜線を通る

選択肢が増えます.
ある点に行く点は最大この3つしかないのですから, この3つのうち最も小さいcostが, この点にたどり着く最小costと考えられます. このとき, 上2つはcost 1がかかるので, +1しますが, ボーナス経路はcost 0なのでそのまま移してきます.

これをすべての点について繰り返すと
f:id:Constellation:20091021161400p:image
以上のようになります.
そして当初の定義どおり, 右下に来たときの最小costが編集距離なのですから, 当然編集距離は4と導かれます.

Diff

これで編集距離およびlevenshtein距離のgraphが求まりました. このとき実はもうすでにdiffも求まっているようなものです.
ということでdiffを求めます. 以下注意点.

  • このとき, 経路を先頭からたどると, 当然途中で別経路のより少ない処理の経路に捨てられ, 行き止まりが多くなります.
  • しかし, よく考えれば, 後ろからたどれば必ず先頭まで到達できるわけです.

後ろから, より少ない編集距離をとる数値の経路を通ります.
f:id:Constellation:20091021161401p:image
では実際にJSで(何で? 好きだからです)実装してみましょう.
肝は斜線の考え方ですが, あらかじめ引いてあると考えるより, 文字列が同じときだけ候補が増えるという風にすると簡単になります. コメントをつけているのでそこも参照してみてください.

// edit graph
function diff1(a, b){
  var a_list = a.split('');
  var a_len = a_list.length;
  var b_list = b.split('');
  var b_len = b_list.length;

  // levenshtein
  var grid = new Array(a_len+1);
  for(var i = 0;i<=a_len;++i){
    grid[i] = new Array(b_len+1);
    for(var j = 0; j <= b_len; ++j){
      grid[i][j] = 0;
    }
  }
  for(i = 1; i <= a_len; ++i){
    grid[i][0] = i;
  }
  for(j = 1; j <= b_len; ++j){
    grid[0][j] = j;
  }
  var row, col;
  for(row = 1; row <= a_len; ++row){
    for(col = 1; col <= b_len; ++col){
      var cost_a = grid[row-1][col]+1;
      var cost_b = grid[row][col-1]+1;
      var min = Math.min(cost_a, cost_b);
      if(a_list[row-1] === b_list[col-1]){
        var cost_c = grid[row-1][col-1];
        min = Math.min(min, cost_c);
      }
      grid[row][col] = min;
    }
  }
  // grid => levenshtein result

  // snake
  var result = [];
  row = a_len;
  col = b_len;
  // 後ろから見ると, 行き止まりなしでたどれる
  while((row > 0) && (col > 0)){
    // ある位置にいけるのは上, 隣, 斜め上の3方向
    // Edit graphの定義より, 下への移動はdelete, 横への移動はadd
    // 斜め上への移動はcommon
    var a = grid[row][col-1];
    var d = grid[row-1][col];
    var c = grid[row-1][col-1];

    if(d < a){
      // しかし, 斜め上に行くには文字列が一致している必要がある
      if((a_list[row-1] === b_list[col-1]) && (c < d)){
        result.push('|  '+a_list[row-1]);
        --row;
        --col;
      } else {
        result.push('- '+a_list[row-1]);
        --row;
      }
    } else {
      if((a_list[row-1] === b_list[col-1]) && (c < a)){
         result.push('|  '+a_list[row-1]);
         --row;
         --col;
      } else {
         result.push('+ '+b_list[col-1]);
         --col;
      }
    }
  }

  // rowが0になってしまって, まだcolが残っているとき
  // Edit graph最上段で横にずれ続ける
  while(col > 0){
    result.push('+ '+b_list[col-1]);
    --col;
  }

  // colが0になってしまって, まだrowが残っているとき
  // Edit graph最左端で上へずれ続ける
  while(row > 0){
    result.push('- '+a_list[row-1]);
    --row;
  }

  // show result
  return result.reverse();
}

var result = diff1("string", "strength");

for(var i = 0, len = result.length; i < len; ++i)
  console.log(result[i]);

実行結果

| s
| t
| r
- i
+ e
| n
| g
+ t
+ h

こうして最も単純な考え方でのdiffを求めました.

Myers氏らの手法による(the O(ND) algorithm)によるdiff

概要

上のやり方でdiffをとれました. しかしよく考えると右上や左下など, 到底使わないだろう辺境の距離まで求めてしまっています. これでは時間がかかってしまいます. それを解決するのがMyers氏らの計算量O(ND)のalgorithmです. このalgorithmは方眼紙を用意して自分で書いてみると非常にわかりやすくなります.

Myers氏らのalgorithmの決定的な違いは, 編集距離ごとに考えている点です. 上の場合はとりあえず右下まで行っていましたが, 今度は編集距離(以下D)ごとにいけるかどうかを試しています. つまり,

D = 0 でいけるかなー? => 無理!
D = 1 でいけるかなー? => 無理!
D = 2 で(ry           => 無理!
D = 3 で(ry           => 無理!
D = 4 で(ry           => いけた! => D = 4だぜ!

という風な考えをとります.

D=0とk線の導入

このとき, 新たな考え方を導入します. これが自身が2日悩みまくったk線というものです.
f:id:Constellation:20091021161359p:image

このように斜めにk線を引きます. また, そのk線に番号をつけて識別しています. この番号は
kの番号 = x座標 - y座標
という考えでつけています.{(0, 0)スタート}
こうして準備が整いました. このとき, 以下の点に注意すると理解しやすいかと思います.

  • 編集距離最短とは, つまりcost 0である斜線をいかに通りまくるかということにかかっている.

D=0

まず D = 0のときを考えましょう.
D = 0 のとき, 縦にも横にも移動するにはcost 1がかかるので動けません. つまり
「斜線だけでいけ!!!!」
ということになります. で斜線だけで最高でどこまでいけるかを考えるわけです.
今回はstrまで共通なため,
f:id:Constellation:20091021161403p:image
というふうに(0, 0)から(3, 3)まで行くことができます. しかしこれ以上動けません. よって D = 0 では到達できないということがわかります.
このとき, k = 0 の線を通りましたね. そこで, k = 0の線では(3,3)までいけたということを記録しておきます. このとき, メモリにうるさいalgorithmはこういいます.

kの番号 = x座標 - y座標
なのだから k = 0のときに y = 3までいけたとしておけばxは記録しなくていい...

D=1

さあ, 次にD = 1のときを考えます. 一回だけ縦か横に動いていいわけです.
さてここからが肝.
先ほどk = 0のときに y = 3までいけたのですから(これを思い出すために記録しておかないといけない)そこからはじめるのが効率がいいでしょう.
このときcost 1だけ動くと, 対象はk = -1の線かk = 1の線上に乗ることになります.
f:id:Constellation:20091021161404p:image
つまり, D = 1のとき, 利用できるk線は -1 <= k <= 1 なのです.
これを一般化すると

D = n のとき -n <= k <= nのk線が利用できる

ということになります.
かつ, 気をつけるべきこととして, 枠外に出てしまっては意味がありません. よって, k線の範囲は
(-文字列Bの長さ) <= k <= (文字列Aの長さ)
である必要もあります. この場合は, "string"で6, "strength"で8より

 -8 <= k <= 6

です. 図を見てもここまでしか考慮されていないことがわかります.

次に, さらに逆に考えて見ます. k線ごとに考えるわけです.

今 k = -1の線について考えます. D = 1のとき, k = -1の線上に移動できるのは基本的に,
f:id:Constellation:20091021161402p:image

  • k = 0 (現在のk+1の線) のときに行ききった場所から下に1ずれてくる
  • k = -2 (現在のk-1の線) のときに行ききった場所から横に1ずれてくる

の2つの道が考えられます.(右からずれてくるとかは, 右下に近づこうとしている以上ありえません) しかし, k = -1 は現在利用できるk線のなかで端っこであり, 横からずれてくるルートは存在しません. よって必然的に上からずれてくるルートになります.

このように, 条件が限られる場合があります.

  • k = -d, k = -N のとき, 横には存在しないので, 上からのみ移動してくる.
  • k = d , k = M のとき, 上には存在しないので, 横からのみ移動してくる.

よってk = -1のときのスタート点はk = 0の点(3,3)からyがひとつずれて(3,4)です.
この(3,4)からD = 0のときと同じように斜めにいけるだけいきます. まあ今回は斜めにまったくいけないので, k = -1のときの最遠端は(3,4)となります. これをまた記録しておきます.

次にD = 1のときのk = 0について考えるといきたいところですが, その必要はありません. なぜならD=0で利用できる点から1移動して, k = 0線上に乗っかる候補はありえないからです. つまり, D = 1のときには k = -1 と k = 1のみ調べればいいわけです.

これを一般化すると

D = n のとき, -n <= k <= nのk線が利用でき,
また, それは -n, (-n+2*1), (-n+2*2), ... , n
という風に2ずつ増えていく.

という風にすることができます. ソースコード的にいうと

for(var k = (-D); k <= D; k+=2){
  // k線ごととの処理
}

という風にすればいいのが見えてきます.

では同じようにk = 1について考えて見ましょう.
f:id:Constellation:20091021161406p:image
このとき, k = -1のときに示した2つの経路のうち, 横からずれる経路しかつかえません.
よってk = 1のスタート点はk = 0の点(3,3)からxがひとつずれて(4,3)になります.
当然これももういけないので, k = 1のときの最遠端は(4,3)であると記録しておきます.

D = 1について調べきりましたね. このときも目標点(6,8)にはたどりつきませんでした.

D=2

D = 2のときの調べるべきk線は(-2, 0, 2)ですね. 端っこの-2, 2についてはk = 1, -1と同じようにすればいいので, k = 0のときを考えます.
k = 0のとき,
f:id:Constellation:20091021161407p:image

  • k = 1 (現在のk+1の線) のときに行ききった場所から下に1ずれてくる
  • k = -1 (現在のk-1の線) のときに行ききった場所から横に1ずれてくる

の2つとも考えることができます. ここでどちらをとるべきか考えましょう.
横にずれてくるとき, 至る座標は(3,4)にxを1足して(4,4)です.
下にずれてくるとき, 至る座標は(4,3)にyを1足して(4,4)です.
今回は同じなのでどちらをとってもいいでしょう.
しかし, これがずれているとき, 目標としては(6,8)に近づくのですから, それに近いほう, つまりより数値の大きいほうをとればいいのです. 同じcostでより進んだほうをとるわけです.

さあどこまでいけるでしょうか. 斜線部分を通ることができるので, 一気に(6,6)まで進むことができました. k = 0のときはy=6, これを記録します.
下のグラフに, k = -2, k = 2のときのいける距離も記述しておきました.
f:id:Constellation:20091021161405p:image
結局D=2ではまだ(6,8)には到達できていません. しかしだいぶ近くなりました.

D=3

D = 3のときの調べるべきk線は(-3, -1, 1, 3)ですね. これを先ほどと同じようにやると以下のグラフのようになるのがわかるでしょうか.
f:id:Constellation:20091021161409p:image
ひとつ注意してほしいのは(6,6)までいった k = 0 の線からの移動です.
もう端まで来ているので, 横には動けなくなっていると考えるかもしれません. k = 1 のとき, 横からスライドする(7,6)でなく, 上から降りてくるルートを取るようにするのが自然です. しかし, 実は隣から行くルートをとっても問題ありません.
なぜなら, 枠外には斜線が存在しないので, これからこのルートが, 最短になるということはまずなくなってしまうからです. つまり, 計算しても意味がないことが確定している部分となります.

これでも到達できていません. しかし後一歩です. もうすでに処理が反復になっています. 新しいことは出てきていません.

D=4

D = 4のときの調べるべきk線は(-4, -2, 0, 2, 4)ですね.
このとき, k = -2の線について注意してください. k = -1 or k = -3から動いてくるわけですが, 最も(6,8)に近い(6,7)からの移動によって, ついに(6,8)に至るわけです!!!!
f:id:Constellation:20091021161410p:image
やりました! こうして経路および D = 4の値が確定しました!!!!!!!!!!!!!!

Diff

さあDiffを出しましょう. Levenshteinのときに, 後ろからだと行き止まりなしだったのは, この場合も成立します.
つまり
f:id:Constellation:20091021161408p:image
のような木構造を作っていけばいいわけです. そして最後にたどり着いたときの木の最後尾からたどっていけば, ルートが取れます.
それらすべてを考慮したJSのソースを以下に記述します.

// Myers's algorithm
function diff2(str_a, str_b){
  var A = str_a.split('');
  var B = str_b.split('');
  // str_aの長さM
  var M = A.length;
  // str_bの長さN
  var N = B.length;

  // それぞれのk線上で最も進んだ結果の位置
  // およびそこにいたるまでの経緯が保存される配列
  var V = [];

  // offsetは, 下で使うkが最低で-Nになることから,
  // offsetを足すことで配列のindexの先頭, 0にする目的がある.
  // kの範囲は最低-Nのため, k+offset=0となるよう設定
  var offset = N;

  // Vの初期化
  // M + N + 1の1 は, 原点を通っている線が数えられているから
  for(var i = 0, len = M + N + 1; i < len; ++i){
    V[i] = {
      // y座標.
      // kが与えられることから, xはkとyとでもとめられるので覚える必要なし
      y    : 0,
      // tree先頭 ( というか最後からたどるので, 最後 )
      // これをたどっていくことで, 後に編集履歴を参照できる
      tree : null
    }
  }

  found: {
    // currentは現在のk線上の結果として保存すべきVの値
    // 周りを調べて, この与えられたcurrentに適切な値を入れて回る
    var current = null, add = null, del = null;
    for(var d = 0, len = M + N + 1; d < len; ++d){
      // ひとつずつDがずれると, kを2つずつ移動することで
      // kに関する部分のみ考えることができる
      for(var k = -d; k <= d; k+=2){
        // 枠外に出てしまっている
        // k の範囲は最大 -N <= k <= M
        if((k < -N) || (M < k))
          continue;

        current = V[k+offset];

        // 一番初め, d === 0のときは, どちらからも進んでくることができない
        if(d !== 0){
          if(k === -d || k === -N){
            // 隣からくるルートは存在しないので,
            // 自然と上から降りてくるルート
            // 降りてくるのでy座標は1増える
            add = V[k+1+offset];
            current.y = add.y + 1;
            current.tree = {
              type : '+',
              prev : add.tree
            }
          } else if (k === d || k === M){
            // 上から降りてくるルートは存在しないので,
            // 自然と隣からスライドするルート
            del = V[k-1+offset];
            current.y = del.y;
            current.tree = {
              type : '-',
              prev : del.tree
            }
          } else {
            // 隣から降りるルートと上から降りてくるルートの二つのうち
            // より下に降りてこられるルートを取る
            if(V[k-1+offset].y > (V[k+1+offset].y+1)){
              // 隣からスライドするルート
              del = V[k-1+offset];
              current.y = del.y;
              current.tree = {
                type : '-',
                prev : del.tree
              }
            } else {
              // 上から降りてくるルート
              // y座標は1増える
              add = V[k+1+offset];
              current.y = add.y + 1;
              current.tree = {
                type : '+',
                prev : add.tree
              }
            }
          }
        }
        // 定義(k = x - y)より
        var x = current.y + k;
        var y = current.y;

        // snake
        while(x < M && y < N && (A[x] === B[y])){
          current.tree = {
            type : '|',
            prev : current.tree
          }
          ++x;
          ++y;
        }
        current.y = y;
        if(x >= M && y >= N){
          // tree 完成
          // 一応編集距離dをセットしておいた
          current.d = d;
          break found;
        }
      }
    }
  }

  // tree 完成後処理
  // treeの終端から先頭に向かってたどって行く
  // prevをたどって後ろから戻して行く
  var list = [], type = null, a_index = M - 1, b_index = N - 1;
  for(var i = current.tree; i != null; i = i.prev){
    // "+" or   "-"  or "|"
    // add or delete or common
    type = i.type;
    if(type === '+'){
      // addのとき, Bを調べて足された文字を調べる
      list.unshift(type+' '+B[b_index]);
      --b_index;
    } else if(type === '-'){
      // delのとき, Aを調べて削られた文字を調べる
      list.unshift(type+' '+A[a_index]);
      --a_index;
    } else {
      // commonのとき, どちらも同じなのでとりあえずAを調べている
      list.unshift(type+' '+A[a_index]);
      --a_index;
      --b_index;
    }
  }
  return list;
}
var result = diff2("string", "strength");

for(var i = 0, len = result.length; i < len; ++i)
  console.log(result[i]);

k線ごとの最も近づいたときのy座標, そして, そこに至るまでの経路木を保存し, 木を後ろから付け足しながら動いているのがわかるでしょうか.
JSはGCが走るのでこんな宙にういたような木構造によってDiffをとっていますが, C, C++などではそうはいきません. その際は, 配列にNodeをつめ続け, k線ごとの値には, その配列のindexを保存しておくという考え方がきれいかもしれません.
C/C++については, わかりやすくすばらしいものがあるので, ぜひそちらを参照してください.
2009-04-27 - エンジニアのソフトウェア的愛情

Wu氏らの手法(the O(NP) algorithm)によるdiff

概要

このalgorithmは先ほどのものとは違い, 今度は損失を考えるやり方です.
まずはじめに, 長さM文字の文字列Aと長さN文字の文字列Bでは最低でも
編集距離D = |M - N|
がかかってしまいます. つまり, D < |M - N|について調べていたのは本質的には無駄であったということです. しかし, 先のalgorithmはDが一つ少なかったときの数値が必要になり, どうしても計算しないといけません.
しかしこれを解決する方法があります. それがWu氏らのalgorithmです.

最小編集距離 delta

k = x - yと前のように定義し, k 線を引きます.
f:id:Constellation:20091021162812p:image
このとき, 変換元文字列Aの長さはBよりも長いものであるという前提が必要となります. これは, 長さ M, Nを比較し, 場合によってはひっくり返して, すべての処理を逆に記録すればいいだけなので問題ありません.

delta = M - N
とします. M > N になるようにひっくり返し済みなので, 必ず正になります. このdeltaが編集距離としてありうる最小の値です.
このとき,
f:id:Constellation:20091021162813p:image
の部分が 編集距離D(以後Dで統一) = delta で到達可能であるときに通る可能性のある経路部分です. 0 <= k <= delta の部分ですね.
さあ本題に近づくわけですが,
最終目的である(8, 6)は, k = delta 線上に存在している.
ということを頭に入れておくと考えやすいかもしれません.

損失p

Myers氏らのalgorithmのときは, 編集距離Dに着目していましたが, 今度は損失pについて着目します.
述べたように, 最もうまくいけばD = deltaで編集を終えることができますが, 実際には斜線が足りないのでいけないわけです. そこで, 本来deltaでいけるところを迂回してしまったと考え, その迂回する回数を, deltaからの損失pとして考えます.
考えやすいように小さなものを下に示します.
f:id:Constellation:20091021163834p:image
このように, 色のついた部分から逸脱し, 迂回していますね. このとき距離1迂回しているので p = 1とします. これでたどり着いたとすると,
D = (M - N) + p * 2
つまり, D = (3 - 2) + 1 * 2 = 3 となります.
なぜなら, 迂回すると帰って来るためのcostも必要となるため, p * 2のcostが余分にかかるわけです.

この損失pについて考えます. つまり, "string" と "strength"のとき,

p = 0 でいけるかなー? => 無理!
p = 1 でいけるかなー? => いけた! => D = (8 - 6) + 1 * 2 = 4だぜ!

という風な考えをとります.

迂回の考え方

上で, さらっと迂回するといいましたが, 迂回の考え方は場所によって異なります.
f:id:Constellation:20091021163557p:image
2とおりの迂回ルートが描かれていますが, 迂回第一歩目が横と縦で違うのがわかるでしょうか. 一般に, k線が

  • -p <= k < delta のとき
  • k = delta のとき
  • delta < k <= delta + p のとき

の3つに分類することができます.

-p <= k < delta のとき, はじめの点を取る経路は2つ存在します.

  • k-1の線から横にひとつずれてくる
  • k+1の線から縦にひとつずれてくる

k - 1の線からずれるのは, k = delta に近づく経路なので, pは現在と同一です.
しかし, k + 1の線からずれてくるのは, k = delta から遠ざかる経路なので, p はk+1のとき記録した値より+1されたものになります. つまり, 現在のpにあわせるには, 一つ少ないp-1のときに記録した経路をとってこないといけないわけです.
つまり,
初回位置の点のy座標をS,とするとき,

S = max(pが同じときのk-1のk線上のy座標+1, pが1少ないときのk+1のk線上のy座標)

となります.

k = delta のとき, はじめの点を取る経路は2つ存在します.

  • k-1の線から横にひとつずれてくる
  • k+1の線から縦にひとつずれてくる

k - 1の線からずれるのは, k = delta に近づく経路なので, pは現在と同一です.
k + 1の線からずれるのは, k = delta に近づく経路なので, pは現在と同一です.
つまり,
初回位置の点のy座標をS,とするとき,

S = max(pが同じときのk-1のk線上のy座標+1, pが同じときのk+1のk線上のy座標)

となります.

delta < k <= delta+p のとき, はじめの点を取る経路は2つ存在します.

  • k-1の線から横にひとつずれてくる
  • k+1の線から縦にひとつずれてくる

k - 1の線からずれるのは, k = delta から遠ざかる経路なので, pは現在-1です.
k + 1の線からずれるのは, k = delta に近づく経路なので, pは現在と同一です.
つまり,
初回位置の点のy座標をS,とするとき,

S = max(pが1少ないときのk-1のk線上のy座標+1, pが同じときのk+1のk線上のy座標)

となります.

何いってんだろうってなるのではないでしょうか. 自分だったら「はあ?」ってなります.

pの迂回解説

意味わかんない的状態になりかねないので書きまくります. 上の説明でわかった方は飛ばしてもらえれば.

今, 図のようなEdit Graphを考えます.
f:id:Constellation:20091021163558p:image
M = 6, N = 4, delta = 2
最小編集距離として期待できるのは D = 2ですね.

これについてp = 0で考えます.
k=0のとき, 始点Sは(0,0)です. ここからまず, いけるだけ斜めに行きます.
f:id:Constellation:20091021163556p:image
よって, 図のように(2,2)まで行きます. ここで,
k = 0, p = 0のとき(2,2)
と記録します. この(2,2)は何を表しているのでしょうか.

ここからどううまいこといったとしても, 目的地(6,4)にいたるにはあとcost 2必要になります(横に2cost).
それを現在期待しているD = 2から引いた0がこのk=0において使うことのできる編集距離であるということです. 今回編集距離0によってk=0のとき(2,2)までこれました. この(2,2)はつまり, D=2のときk=0において来ることのできる最遠端であり, p=0のときにk=0において来ることのできる最遠端でもあります.

次にk=1のときを考えます. k = 1は -p <= k < delta なので, 次の2つ

  • k-1の線から横にひとつずれてくる
  • k+1の線から縦にひとつずれてくる

の経路を通って始点Sにいたります.
このとき, 前に述べたように, k+1の線からずれてくるときは, p = -1(現在のp-1)のときの記録でなければなりません. これは何故でしょうか.

先ほどk=0のときp=0において来ることのできる最遠端を求めていました. これは, 残りの経路がすべて目的地に近づくように移動する. つまり残っているD = 2はすべて近づくために横移動すると考えられたうえでの値です. しかし, このとき, 次に下に移動するとどうでしょうか.
目的地のk線, k = deltaから離れていますね. つまり迂回しているわけです. よってその点に行くには, p = 0では行くことができません. その点においては迂回してしまったので p = 1となっているのです.

これを逆から見ると理解できます. k = 1のときの話に戻って, S に上から降りてくるルートを取るとき, 迂回してその点まで来ているわけです. よってSに来るときに損失がp であるためには, 迂回前の点では損失がp-1でなければいけないのです.

このことを考えると, 上の迂回の際の理由がわかるかもしれません. 説明下手ですいません... 下のループと実例見たほうが速いかもしれません.

ループの順番

始点Sを求めるために異なるpの2点が必要になっています. では, k と p の2次元配列にするのでしょうか?
しかし, うまくすれば1次元でいけます.

まず, -p <= k < delta のとき,

S = max(pが同じときのk-1のk線上のy座標+1, pが1少ないときのk+1のk線上のy座標)

でしたね. このため, 配列の更新を-pからdeltaの順番ではじめていけば, k = nの値を更新するとき, k=n+1はまだ更新していないのでp-1のときの値, k=n-1の値は更新済みなのでpのときの値になっており, ちょうど良くなります.
このため, はじめのループは

for(k = -p; k < delta; ++k){
  // ...
}

とすればいいと考えられます.

次に, k = deltaのとき,

S = max(pが同じときのk-1のk線上のy座標+1, pが同じときのk+1のk線上のy座標)

でしたね. どちらもpについてなので, 更新済みの情報が必要です. よって, k = deltaは他のkがすべて更新が終わってからしか調べることができないのがわかるでしょうか.

そして, delta < k <= delta+p のとき,

S = max(pが1少ないときのk-1のk線上のy座標+1, pが同じときのk+1のk線上のy座標)

でしたね. このため, 配列の更新をdelta+pからdeltaまで, -1ずつ下げながらすれば, k=nの値を更新するとき, k=n+1は更新済みなのでpのときの値になっており, k=n-1の値はまだ更新していないのでp-1のときの値になっていてちょうど良くなります.
このため次のループは

for(k = delta+p; k > delta; --k){
  // ...
}

とすればいいと考えられます.

以上をまとめれば, 次の順番でループをまわせば, 1次元配列で処理を終えることができるとわかります.

// -p <= k < delta
for(k = -p; k < delta; ++k){
  // ...
}

// delta < k <= delta+p
for(k = delta+p; k > delta; --k){
  // ...
}

// k = delta
k = delta;
// ...
demo

"string" と "strength"の例でやってみます.
まず, "string" は長さ6, "strength"は長さ8より, この対象を入れ替えます.
つまり,

文字列A "strength" M = 8
文字列B "string"   N = 6
delta = 8 - 6 = 2
入れ替えられてる(exchanged = true)

という条件です.
この条件でEdit Graphに表すと, 以前のGraph
f:id:Constellation:20091021162813p:image
となります. さあstartです.

p=0

p=0のときを考えます.
f:id:Constellation:20091021163827p:image
k=0のとき, 始点Sは(0, 0)となります. 始点を求めようにも両側とも未定義ですので. そして斜めに行けるだけ行くと, (3,3)までいくことができます. k = 0, p = 0のとき(3,3)ですね. これをfp[0,0] = 3(y座標) と表しておきます.

先ほどのループの順番に乗れば次はk=1です.
f:id:Constellation:20091021163828p:image
k=1のとき, 始点Sはk-1の線のp=0のときの値, つまり先ほどの(3,3)か, k+1, p=-1のときの値から求めるわけですが, k=2, p=-1の値fp[2,-1]は未定義ですので自動的に(3,3)から横にスライド, S=(4,3)です. ここから斜めには進めないので, fp[1,0] = 3です.(今のところ, fpを2次元配列のように書いているが, 上で述べた説明のとおり, 1次元で表すことができる)

ここまでで最初のループは終了ですね. 次はdelta+pからdeltaより上まで下がりながら見るわけですが, 現在p=0より, 2番目のループは何もしません.

なので, 次はk=deltaの番です.
f:id:Constellation:20091021163826p:image
k=2のとき, 始点Sはk-1の線のp=0のときの値, つまり前の(4,3)か, k+1, p=-1の...という風に前と同じくなりますが, p=-1のときはまだ未定義なので, 自動的に(4,3)から横にスライド, S=(5,3)です. またこれも斜めにはいけないので, fp[2,0] = 3です.

この時点でfp[delta, p]の値はfp[2,0]より3です. これがN, つまり6になったら目的地に到達ですが, 今のところはまだです.

p=1

p=1のときを考えます.
ループの順番的に, -1..k..1 をみて, 3..k..3 をみて, k = 2を見ればいいみたいですね.
f:id:Constellation:20091021163830p:image
k=-1のとき, 始点Sを, k=-2, p=1から(fp[-2,1])か, k=0, p=0(fp[0,0])からか求める必要があります. fp[-2,1]は未定義なので, fp[0,0]の値(3,3)から縦に降りたS=(3,4)です. ここから斜めにはいけないので, fp[-1, 1] = (3,4)です.
f:id:Constellation:20091021163831p:image
k=0のとき, 始点Sを, fp[-1,1]からか, fp[1,0]からか求めます. 今回初めてどちらも使えますね. fp[-1,1]は先ほどより(3,4), fp[1,0]は前に求めた(4,3)です. fp[1,0]は下に下りてくると(4,4)になります. この場合どちらのyも4なのでどちらでもいいですが, delを先にとるように後で実装したので, とりあえずdel方向, つまり横スライドを取りましょう. S=(4,4)です. ここから斜めに進めるだけ進むと(6,6)にいけます. よってfp[0,1] = (6,6)です. かなり来ましたね.
f:id:Constellation:20091021163832p:image
k=1のとき, 始点Sを, fp[0,1]からか, fp[2,0]から求めます. fp[0,1] = (6,6), スライドして(7,6) fp[2,0] = (5,3), fp[2,0]はひとつ降りて(5,4), y座標を比較すれば(7,6)をとりますね. S=(7,6)です. もう斜めにはいけないのでそのまま. fp[1,1] = (7,6)です.

ここでk=2直前まで来ました. k=delta直前です. 次はk=3ですね.
f:id:Constellation:20091021163833p:image
k=3のとき, 始点Sを, fp[2,0]からか, fp[4,1]からかも止めます. k=2からk=3への移動はk=deltaから遠ざかるので迂回よりpが一つ少ないp=0, k=4からk=3への移動は近づくのでp=1のままでよし. fp[4,1]は未定義なので, fp[2,0]=(5,3)から横にスライドしS=(6,3). 斜めにはいけないので, fp[3,1]=(6,3)です.

つぎはk=delta, k=2ですね.
f:id:Constellation:20091021163829p:image
k=2のとき, 始点Sを, fp[1,1]か, fp[3,1]からもとめます. どちらもk=deltaに近づくのでp=1のままですね. fp[1,1]=(7,6)で横スライド(8,6), fp[3,1]=(6,3)で縦スライド(6,4). 選ぶのはyの大きいS=(8,6)ですね. ここから斜めにはいけないので, fp[2,1] = (8,6)です.

さあループ終了につきfp[delta, p], つまりfp[2,1]を見ると(8,6)となんと到達しました! よって編集距離Dは

D = delta + p * 2 = 2 + 1 * 2 = 4

となり, これまでの経路を保存していればきちんとたどることもできます!!!!!!

Diff

経路の保存はMyers氏らのalgorithmと変わりません.

// Wu's algorithm
function diff3(str_a, str_b){
  var A = str_a.split('');
  var B = str_b.split('');
  var M = A.length;
  var N = B.length;

  var exchanged = false;

  // 長さによって入れ替え
  // ひっくり返す.
  if(N > M){
    exchanged = true;
    var t;
    t = A;
    A = B;
    B = t;

    t = M;
    M = N;
    N = t;
  }

  var offset = N;
  var delta = M - N;
  var size = M + N + 1;
  for(var i = 0, fp = []; i < size; ++i){
    fp[i] = {
      // y座標.
      // kが与えられることから, xはkとyとでもとめられるので覚える必要なし
      y    : null,
      // tree先頭 ( というか最後からたどるので, 最後 )
      // これをたどっていくことで, 後に編集履歴を参照できる
      tree : null
    }
  }
  // snake処理
  // 斜線部分をいけるだけ行く
  var snake = function(k){
    // if分岐は, わざと冗長な書き方をしています
    if((k < -N) || (M < k)){
      return;
    } else {
      var current = fp[k+offset];
      if(k === -N){
        // 横が存在しない
        // down  => 上から降りてくる => ADD
        // down の y座標は上から降りてくるので+1
        var down  = fp[k+1+offset];
        current.y = down.y+1;
        current.tree = {
          type : '+',
          prev : down.tree
        }
      } else if(k === M){
        // 上が存在しない
        // slide => 横からスライド   => DEL
        var slide = fp[k-1+offset];
        current.y = slide.y;
        current.tree = {
          type : '-',
          prev : slide.tree
        }
      } else {
        var slide = fp[k-1+offset];
        var down  = fp[k+1+offset];
        if(slide.y === null && down.y === null){
          // どちらも未定義 => (0, 0)について
          current.y = 0;
        } else if(down.y === null || slide.y === null){
          // どちらかが未定義状態
          if(down.y === null){
            current.y = slide.y;
            current.tree = {
              type : '-',
              prev : slide.tree
            }
          } else {
            current.y = down.y+1;
            current.tree = {
              type : '+',
              prev : down.tree
            }
          }
        } else {
          // どちらも定義済み
          if(slide.y > (down.y+1)){
            current.y = slide.y;
            current.tree = {
              type : '-',
              prev : slide.tree
            }
          } else {
            current.y = down.y+1;
            current.tree = {
              type : '+',
              prev : down.tree
            }
          }
        }
      }
      var y = current.y;
      var x = y + k;
      while(x < M && y < N && A[x] === B[y]){
        current.tree = {
          type : '|',
          prev : current.tree
        }
        x++;
        y++;
      }
      current.y = y;
    }
  }

  var p = -1,
      k =  0;
  // k = deltaのときのy座標が fp[delta+offset]
  // 目的地はk = delta上にあるので,
  // fp[delta+offset]のy座標がNであるとき, x座標もMであり, 目的地
  while(fp[delta+offset].y < N){
    p = p + 1;
    for(k = -p; k < delta; ++k){
      snake(k);
    }
    for(k = delta+p; k > delta; --k){
      snake(k);
    }

    snake(delta);
  }

  var current = fp[delta+offset];
  current.d = delta + 2 * p;
  // tree 完成後処理
  // treeの終端から先頭に向かってたどって行く
  // prevをたどって後ろから戻して行く
  var list = [], type = null, a_index = M - 1, b_index = N - 1;
  for(var i = current.tree; i != null; i = i.prev){
    // "+" or   "-"  or "|"
    // add or delete or common
    type = i.type;
    // exchangedも考慮するように
    if(type === '+'){
      // addのとき, Bを調べて足された文字を調べる
      // ひっくり返されてたら, delにする
      list.unshift(((exchanged)? '-' : '+')+' '+B[b_index]);
      --b_index;
    } else if(type === '-'){
      // delのとき, Aを調べて削られた文字を調べる
      // ひっくり返されてたら, addにする
      list.unshift(((exchanged)? '+' : '-')+' '+A[a_index]);
      --a_index;
    } else {
      // commonのとき, どちらも同じなのでとりあえずAを調べている
      list.unshift('| '+A[a_index]);
      --a_index;
      --b_index;
    }
  }
  return list;
}

var result = diff3('strength', 'string');
// 結果表示
for(var i = 0, len = result.length; i < len; ++i)
  console.log(result[i]);

終了

以上, 正攻法の出し方, Myers氏らの計算量O(ND)のalgorithm, Wu氏らの計算量O(NP)のalgorithmの編集距離導出および, Diffのとり方となります.
むずかしい... というか考えた人すごすぎる.

参考文献

Myers氏らの論文
An O(ND) Difference Algorithm and Its Variations (PDF)
Google Doc
中に, 擬似コードがそのまま載ってたりしたりの充実の内容

Wu氏らの論文
An O(NP) Sequence Comparison Algorith (PDF)
Google Doc

viviの作者による有名どころ
•¶‘”äŠridiffjƒAƒ‹ƒSƒŠƒYƒ€

/0 » diff(1) [Edit Graph]
/0 » diff(2) [Myers' algorithm]
/0 » diff(3) [Wu's algorithm]

C++ Diff Source
2009-04-25 - エンジニアのソフトウェア的愛情 [Edit Graph Diff]
2009-04-27 - エンジニアのソフトウェア的愛情 [Myers' algorithm]

JS Wu's algorithm Diff
diff O(np) javascript implementation « ku

補足

間違いの指摘などありましたら是非是非お願いします.