Markov Cluster Algorithm を Ruby & NArray で実装

NGSデータを用いたde novo assemblyにはクラスタリングアルゴリズムが重要であろう，という妄想からふらふらとクラスタリング法を調べていて見つけたのがMCLことMarkov Cluster Algorithm．2000年に発表された手法らしく特に新しい分けではないが，それなりに性能は良いらしいとか．

詳しい説明は本家ホームページ（MCL - a cluster algorithm for graphs）に丸投げしさっそく実装を行ってみた．こちらのサイトに（2010-06-28 - (setf yaruki nil) - nlpyutoriグループ）rubyの実装例があったのだが，NArrayを使いたい！と思いNArrayバージョンとして実装してみた．

タブ区切りで，edgelistを作り第一引数に，全ノード数を第二引数に与えて実行するというダサい形．最終的にはm1行列からクラスター情報を取り出す必要があるが，まだ実装しておりません．．下記テストデータを与えると瞬間でクラスタリングが終了！多分まともな答えが出ている気がします．

#!/usr/local/ruby1.9.2/bin/ruby19
require "rubygems"
require "narray"

def main
  # parameters
  file = ARGV.shift
  size = ARGV.shift.to_i
  r    = 2.0
  th   = 0.00001
  
  # read graph
  m1 = NMatrix.float(size, size)
  File.foreach(file) do |line|
    i, j = line.chomp.split("\t")
    m1[i.to_i, j.to_i] = 1
    m1[j.to_i, i.to_i] = 1
  end
  m1.I
  m1.div!(m1.sum(0))
  
  # MCL
  loop do
    m2 = m1**2                              # expansion
    m1 = NMatrix.refer(NArray.refer(m2)**r) # inflation
    m1.div!(m1.sum(0))                      # inflation
    dif = (NArray.refer(m1-m2)**2).sum
    break if dif < th
  end
  
  # output
  p m1
end

main if __FILE__ == $0

テストデータ．こちらを参考にさせていただきました（Fig. 2 http://www.matlab.nitech.ac.jp/~matsuo/NL05-1.pdf）．