Rのanimationパッケージとigraphパッケージでプリム法をアニメーション表示

せっかくさっき(Problem 107 - もうカツ丼でいいよな)プリム法を実装したので,Rのanimationパッケージでクラスカル法のシミュレーション - ぬいぐるみライフ?の真似してアニメーション化してみた.プリム法についてはプリム法(最小全域木問題)を参考にしました.

例示したグラフはProblem 107 - PukiWikiの例題のものを使用.同じ結果になってるのが分かると思う.
プリム法はそれまでに確定しているノードのいずれかに最小のコストで到達できるノードから順に確定していく.確定したノードへの最短到達コストはノードが新たに確定する度に計算し直される.その様子が見て取れる.
以下ソース.色々手間取って思った以上に長くなった.

library(animation)
library(igraph)

## 新たに確定したノードから直接到達可能なノードへのコストを更新
update.m.d <- function(min.dist, dec.velt, graph, dec.flag){
  comp <- min.dist > graph[dec.velt, ] # 新たに確定したノードから各ノードへの距離と比較
  ## 新たに確定したノードからの到達距離の方が短い部分を置き換える(未確定のノード)
  min.dist[comp & !dec.flag] <- graph[dec.velt,][comp & !dec.flag]
  return(min.dist)
}

## Prim's algorithm
prim <- function(d.mat){                # グラフは距離行列として与える
  ## 接続していない部分を0に置き換えて書き直し
  d.mat[d.mat=="-"] <- 0
  d.mat <- matrix(as.numeric(d.mat), nrow=sqrt(length(d.mat)))
  ## 距離行列からグラフオブジェクト作成
  g <- graph.adjacency(d.mat, mode="undirected", weighted=TRUE)
  ## ノードのレイアウト(プロット時に使用)
  lay <- layout.fruchterman.reingold(g)
  ## ノード数
  n <- sqrt(length(d.mat))
  ## color, size
  V(g)$size <- 25                       # ノードのサイズ
  V(g)$color <- "lightblue"             # ノードの色
  V.name <- letters[1:n]                # ノードのラベル
  E(g)$color <- "gray"                  # エッジの色
  E(g)$width <- 1                       # エッジの太さ
  par(cex=1.5)                          # 文字サイズ

  ## 初期状態を出力
  plot(g, layout=lay, vertex.label=V.name, edge.label=E(g)$weight)

  ## Prim法 & 経過出力
  d.mat[d.mat==0] <- Inf           # 接続していない部分をInfに書き直し
  min.dist <- rep(Inf, sqrt(length(d.mat))) # 確定ノードから各ノードへの最小コスト
  dec.flag <- rep(FALSE, length(min.dist))  # 確定フラグ
  min.dist[1] <- 0                          # 初期ノード確定
  dec.velt <- 1                             # 一つ前に確定したノード
  V(g)$color[1] <- "pink"                   # 確定したノードの色を変更
  V.name[1] <- 0                            # 各ノードへの最小コスト
  ## 表示
  plot(g, layout=lay, vertex.label=V.name, edge.label=E(g)$weight)
  mtext("start", side=1, cex=1.5, line=1)
  repeat{                         # 全てのノードが確定するまで繰り返し
    ## 最小コストの更新
    min.dist <- update.m.d(min.dist, dec.velt, d.mat, dec.flag)
    min.dist.l <- length(min.dist[!min.dist==Inf])
    V.name[!min.dist==Inf] <- min.dist[!min.dist==Inf] # 更新結果をノードに表示
    ## 表示
    plot(g, layout=lay, vertex.label=V.name, edge.label=E(g)$weight)
    mtext("update", col="blue", side=1, cex=1.5, line=1)
    ## 確定フラグ
    dec.flag[dec.velt] <- TRUE
    if(sum(!dec.flag) < 1) break        # 全てのノードが確定したら終了
    ## 確定フラグが立っていないノードのうち最も到達コストの小さいもの(の一つ目)を確定する
    dec.velt.old <- dec.velt
    dec.velt <- (1:length(min.dist))[!dec.flag][min.dist[!dec.flag] == min(min.dist[!dec.flag])][1]
    V(g)$color[dec.velt] <- "pink"      # 確定したノードの色を変更
    ## 確定しているノードから今確定したノードへ最短コストでたどり着けるエッジを探す(plot用)
    min.weight <- Inf
    min.e <- numeric(2)
    for(e in (1:n)[dec.flag]-1){        # 確定しているノードのみ
      ## ワンステップで行けないノードはスキップ
      if(length(unlist(get.shortest.paths(g, from=dec.velt-1, to=e))) > 2) next
      ## ノード間の重みが小さければそれを採択
      if(min.weight > E(g, path=c(dec.velt-1, e))$weight){
        min.weight <- E(g, path=c(dec.velt-1, e))$weight
        min.e <- c(dec.velt-1, e)
      }
    }
    E(g, path=min.e)$width <- 3         # 選択されたエッジを太く
    E(g, path=min.e)$color <- "red"     # 同じく色を赤に
    ## 表示
    plot(g, layout=lay, vertex.label=V.name, edge.label=E(g)$weight)
    mtext("select", col="red", side=1, cex=1.5, line=1)
  }
}

## 例示用データ
graph  <- c("-", 16, 12, 21,"-","-","-",
            16,"-","-", 17, 20,"-","-",
            12,"-","-", 28,"-", 31,"-",
            21, 17, 28,"-", 18, 19, 23,
           "-", 20,"-", 18,"-","-", 11,
           "-","-", 31, 19,"-","-", 27,
           "-","-","-", 23, 11, 27,"-")

## gifアニメーションとして出力
saveMovie(prim(graph), interval=2, moviename="prim",
          movietype="gif", outdir=getwd(),
          width=640, height=480)