Rのanimationパッケージでクラスカル法のシミュレーション

Rでクラスカル法を実装してみた．ついでに最小全域木を求める過程をアニメーションで表示させてみた．

クラスカル法についてはこちら．

最小全域木問題(クラスカル法とプリム法)

閉路を作らないように気を付けつつ重みが小さい辺から選んでいく様子が分かる．

ソースは以下の通り．

library(igraph)
library(animation)

# 戻り値はUnion Find用の関数
unite_set <- function(n) {
  p <- rep(-1, n)

  root <- function(x) {
    if (p[x] < 0) {
      return(x)
    }
    else {
      return(p[x] <<- root(p[x]))
    }
  }

  function(x, y) {
    x <- root(x)
    y <- root(y)

    if (x != y) {
      if (x > y) {
        t <- x
        x <- y
        y <- t
      }
      p[x] <<- p[x] + p[y]
      p[y] <<- x
    }

    x != y
  }
}

# クラスカル法
kruskal <- function(v, e) {
  # グラフオブジェクトを作成
  e <- e[order(e[,3]),]
  g <- graph.data.frame(e[,1:2], vertices=v, directed=F)
  lay=layout.fruchterman.reingold(g)

  n <- nrow(v)
  m <- nrow(e)

  # グラフの色やサイズなどを設定
  V(g)$size <- 30
  V(g)$color <- rep("lightblue", n)
  E(g)$color <- rep("grey", m)
  E(g)$width <- rep(1, m)
  E(g)$weight <- e[,3]
  par(cex=1.5)

  # 開始
  plot(g, layout=lay, vertex.label=V(g)$name, edge.label=E(g)$weight)
  mtext("start", side=1, cex=1.5, line=1)

  # Union Findの準備
  u <- unite_set(n)

  # 頂点数 - 1 の数だけ辺を選ぶまでループ
  idx <- 1
  count <- 0
  while (count < n - 1) {
    x <- e[idx, 1]
    y <- e[idx, 2]
    len <- e[idx, 3]

    # 選んだ辺を紫色に
    E(g)$color[idx] <- "purple"
    E(g)$width[idx] <- 3

    # グラフの表示
    plot(g, layout=lay, vertex.label=V(g)$name, edge.label=E(g)$weight)
    mtext(sprintf("check edge: %d - %d: length: %d ", x, y, len), side=1, cex=1.5, line=1)

    selected <- F
    if (selected <- u(x, y)) {
      # 辺を採用
      count <- count + 1
      V(g)$color[x] <- "pink"
      V(g)$color[y] <- "pink"
      E(g)$color[idx] <- "red"
    }
    else {
      # 辺を非採用
      E(g)$color[idx] <- "green"
      E(g)$width[idx] <- 1
    }

    # 辺の採用結果の表示
    plot(g, layout=lay, vertex.label=V(g)$name, edge.label=E(g)$weight)
    mtext(sprintf("check edge: %d - %d: length: %d ", x, y, len), side=1, cex=1.5, line=1)
    mtext(ifelse(selected, "select", "leave"), side=1, line=3, cex=1.5, col=ifelse(selected, "red", "blue"))
    idx <- idx + 1
  }

  # 終了
  plot(g, layout=lay, vertex.label=V(g)$name, edge.label=E(g)$weight)
  mtext("finish", side=1, cex=1.5, line=1)
}

# サンプルのグラフ
v <- data.frame(1:7)
e <- data.frame(matrix(c(
  1, 2, 7,
  1, 4, 5,
  2, 3, 8,
  2, 4, 9,
  2, 5, 7,
  3, 5, 5,
  4, 5, 15,
  4, 6, 6,
  5, 6, 8,
  5, 7, 9,
  6, 7, 11),
  ncol=3, byrow=T))

saveMovie(kruskal(v, e), interval=2, moviename="kruskal",
  movietype="gif", outdir=getwd(),
  width=640, height=480)