- 動的計画法分かんない
- しかも、なんか一次元じゃなくて二次元っぽいから怖い＞＜
と思ってgkbrしていた次回アルゴリズムデザインのゼミで担当の系列アライメント(p246 6.6)のアルゴリズムですが、ミスドでホットカフェオレを3杯くらいおかわりしていたら結構簡単にできました。説明用の資料とかはまだ作ってないけど、たぶんどうにかなるでしょう。ということでとりあえずソース貼っておく。
# -*- coding: utf-8 -*- # アルゴリズムデザイン 6.6 系列アライメント # 動的計画法による最小コストの算出 class SequenceAlignment class Point # バックトラック用のクラス # ポインタがないので。。。 attr_reader :x attr_reader :y def initialize(x = 0, y = 0) @x = x @y = y end end attr_reader :x attr_reader :y attr_reader :delta def initialize(x, y, delta) @x = x.split(//) @y = y.split(//) @delta = delta @a = [] #コストを入れておくための2次元配列 (0..@x.length).each{|i| @a.push [] } end def alignment # initialization (0..@x.length).each{|i| @a[i][0] = i * @delta } (0..@y.length).each{|j| @a[0][j] = j * @delta } # calculate minimum cost for sequence alignment # とりあえず表を埋めていけばいい (1..x.length).each{|i| (1..y.length).each{|j| @a[i][j] = min(alpha(i, j) + @a[i-1][j-1], delta + @a[i-1][j], delta + @a[i][j-1]) } } return @a end def min(a, b, c) if (a < b) && (a < c) return a elsif (b < c) return b else return c end end def delta return @delta end def alpha(i, j) # iとjのindexは1から始まっているので-1しておく a = @x[i-1] b = @y[j-1] vowel = ["a", "b", "c", "d", "e"] if a == b return 0 elsif ( (!vowel.include?(a) && !vowel.include?(b)) || # 母音と母音のマッチング (vowel.include?(a) && vowel.include?(b))) # 子音と子音のマッチング return 1 else # 母音と子音のマッチングコスト return 3 end end def argmin(i, j) a = @a[i-1][j] b = @a[i-1][j-1] c = @a[i][j-1] if (a < b) && (a < c) return Point.new(i-1, j) elsif (b < c) return Point.new(i-1, j-1) else return Point.new(i, j-1) end end def back_track(i = @x.length, j = @y.length, result = []) point = argmin(i, j) if i==0 && j==0 # 最後まで行きついた return result elsif result.length == 0 # 再帰の最初 x = @x[point.x] y = @y[point.y] result.push [x, y] return back_track(point.x, point.y, result) else prev = result[result.length-1] if prev[0] == @x[point.x] x = nil else x = @x[point.x] end if prev[1] == @y[point.y] y = nil else y = @y[point.y] end result.push [x, y] return back_track(point.x, point.y, result) end end def minimum_cost return @a[@x.length][@y.length] end end a = "name" b = "mean" puts "Let us consider sequence alignment betweeen \"#{a}\" and \"#{b}\"." sa = SequenceAlignment.new(a, b, 2) puts "=" * 30 puts "Array for sequence alignment is as follows." sa.alignment.reverse.each{|i| puts "#{i.join(" ")}" } puts "=" * 30 puts "Result is as follows." sa.back_track.transpose.map{|a| a.reverse.map{|b| b.nil? ? "-" : b }.join(" ") }.each{|i| puts i } puts "Minimised cost is #{sa.minimum_cost}."
実行結果。
/Users/syou6162/ruby% ruby sequence_alignment.rb Let us consider sequence alignment betweeen "name" and "mean". ============================== Array for sequence alignment is as follows. 8 6 4 5 6 6 4 4 5 4 4 3 2 3 5 2 1 3 5 6 0 2 4 6 8 ============================== Result is as follows. n a - m e m e a n - Minimised cost is 6.
追記
上のソース、バックトラックの付近がおかしくなってました＞＜。argminのindex、iとjがこんがらかってました。iは縦のほうのindex、jは横のindexになっているので、-1されるものが直感と違ったりしまします。# -*- coding: utf-8 -*- # アルゴリズムデザイン 6.6 系列アライメント # 動的計画法による最小コストの算出 require "pp" class SequenceAlignment class Point # バックトラック用のクラス # ポインタがないので。。。 attr_reader :x attr_reader :y def initialize(x = 0, y = 0) @x = x @y = y end def to_s return "(#{x}, #{y})" end end attr_reader :x attr_reader :y attr_reader :a attr_reader :b attr_reader :delta def initialize(x, y, delta) @x = x.split(//) @y = y.split(//) @delta = delta @a = [] #コストを入れておくための2次元配列 (0..@x.length).each{|i| @a.push [] } @b = [] #バックトラックするための2次元配列 (0..@x.length).each{|i| @b.push [] } end def alignment # initialization (0..@x.length).each{|i| @a[i][0] = i * @delta } (0..@y.length).each{|j| @a[0][j] = j * @delta } (1..@x.length).each{|i| @b[i][0] = Point.new(i-1, 0) } (1..@y.length).each{|j| @b[0][j] = Point.new(0, j-1) } # calculate minimum cost for sequence alignment (1..@x.length).each{|i| (1..@y.length).each{|j| @a[i][j] = [alpha(i, j) + @a[i-1][j-1], delta + @a[i-1][j], delta + @a[i][j-1]].min @b[i][j] = argmin(i, j) } } return @a end def argmin(i, j) # transposeになっていることに注意!! a = delta + @a[i][j-1] # 左 b = alpha(i, j) + @a[i-1][j-1] # 左上 c = delta + @a[i-1][j] # 上 if (a < b) && (a < c) return Point.new(i, j-1) elsif (b < c) return Point.new(i-1, j-1) else return Point.new(i-1, j) end end def delta return @delta end def alpha(i, j) # iとjのindexは1から始まっているので-1しておく a = @x[i-1] b = @y[j-1] vowel = ["a", "b", "c", "d", "e"] if a == b return 0 elsif ( (!vowel.include?(a) && !vowel.include?(b)) || # 母音と母音のマッチング (vowel.include?(a) && vowel.include?(b))) # 子音と子音のマッチング return 1 else # 母音と子音のマッチングコスト return 3 end end def back_track def back_track_prepare(i = @x.length, j = @y.length, result = []) if i==0 && j==0 # 最後まで行きついた return result end point = @b[i][j] result.push point return back_track_prepare(point.x, point.y, result) end prev_x = 1.0/0 prev_y = 1.0/0 str_x = [] str_y = [] back_track_prepare.reverse.each{|i| if i.y == prev_y # 上に移動するとき str_x.push @x[i.x] tmp = str_y.pop # 末尾にあったやつを取り出して str_y.push nil # nilを入れてから str_y.push tmp # 取り出したやつを入れる elsif i.x == prev_x # 右に移動するとき str_y.push @y[i.y] tmp = str_x.pop str_x.push nil str_x.push tmp else # 右上に移動するとき str_x.push @x[i.x] str_y.push @y[i.y] end prev_x = i.x prev_y = i.y } return [str_x, str_y] end def minimum_cost return @a[@x.length][@y.length] end end a = "mean" b = "name" puts "Let us consider sequence alignment betweeen \"#{a}\" and \"#{b}\"." sa = SequenceAlignment.new(a, b, 2) puts "=" * 30 puts "Array for sequence alignment is as follows." sa.alignment.reverse.each{|i| puts "#{i.join(" ")}" } puts "=" * 30 puts "Result is as follows." puts "Minimised cost is #{sa.minimum_cost}." sa.back_track.map{|a| a.map{|b| b.nil? ? "-" : b } }.map{|i| i.join(" ") }.each{|i| puts i }
実行結果。
/Users/syou6162/ruby% ruby sequence_alignment.rb Let us consider sequence alignment betweeen "mean" and "name". ============================== Array for sequence alignment is as follows. 8 6 5 4 6 6 5 3 5 5 4 3 2 4 4 2 1 3 4 6 0 2 4 6 8 ============================== Result is as follows. Minimised cost is 6. m e a n - n - a m e
