このブログの更新は Twitterアカウント @m_hiyama で通知されます。
Follow @m_hiyama

メールでのご連絡は hiyama{at}chimaira{dot}org まで。

はじめてのメールはスパムと判定されることがあります。最初は、信頼されているドメインから差し障りのない文面を送っていただけると、スパムと判定されにくいと思います。

参照用 記事

Rの名前管理とか

Rで離散力学系のアニメーションを作ってみた」で使ったソースコードを“清書”したものをこの記事の末尾に貼り付けます。清書するために、「Rのコーディングスタイルとその背景」を調べたのでした。名前の付け方とコメントの書き方の方針は「Rのコーディングスタイルとその背景」の「当面はこうしようと思う」に書いたとおりです。

Rの場合、いくつかの理由から、「名前の付け方なんて趣味の問題だ」と片づけられないのです。

RのクラスシステムであるS3システムは、ネーミングのコンベンションに依存しています。S3メソッドを書くにはネーミング・コンベンションを守る必要があるし、逆に、S3メソッドではない関数に混同するような名前を付けたくありません。

言語機能として、名前をグループ化して管理する手段がないので、これもネーミング・コンベンションに頼ることになります。パッケージを作るなら名前空間を使えますが、小さなスクリプトファイル程度の規模でパッケージにするのは手間が掛かり過ぎます。

スクリプトファイルという物理的単位は、スコープや名前空間に区切られることがないので、すべてのスクリプトファイルからの名前とライブラリ(パッケージ)からの公開名が単一の大域名前空間に一緒に混ぜてばらまかれます。どの名前がどこから来たか分からず整理整頓するのが大変です。

対話的環境のワークスペース(大域名前空間)が保存されるのはいい点もありますが、ゴミも一緒に保存・再現されるので、R処理系を使い続けるとワークスペースがとっ散らかってドンドン汚くなります。

自前で名前を列挙する関数を付けたのですが、これが良い方法かどうかは分かりません。

#
# 名前の管理
# ==========

# この“モジュール”で使用されている名前を列挙する
ddds.names <- function() {
  list(
    global = c("blur", "blurredPoint", "makeBlurredPointsZMatrix"),
    module = ls(envir = .GlobalEnv, pattern="^ddds\\.*"),
    modulePrivate = ls(envir = .GlobalEnv, pattern="^\\.ddds\\.*", all=TRUE),
    class = ls(envir = .GlobalEnv, pattern="^.+\\.ddds$"),
    classPrivate = ls(envir = .GlobalEnv, pattern="^\\..+\\.ddds$", all=TRUE)
    )
}

名前はすべてキャメルケースを使いました(結局、単語のドット区切りは使わなかった)。関数も変数も小文字始まりですが、ひとつ例外があって、オブジェクトの状態を変更するメソッドの名前は大文字始まりにしました。Reinit.ddds(), DoTrans.ddds() がその例です。

Rのオブジェクトはイミュータブルなので、通常の意味で更新は出来ませんが、新しいオブジェクトを作って返す関数で更新の代用としています。オブジェクトの一部だけを変更するなら、オブジェクトまるまるコピーされるわけではないようなので、効率は悲惨という程じゃないでしょう。

最初は、離散力学系のモデルとなるオブジェクトには、表示やアニメーションのパラメータは入れなかったのですが、結局、dddsオブジェクトに全部のパラメータを詰め込みました。試行錯誤のある時点を保存する機会が多いので、そのときすべてのパラメータが一度に保存されたほうが便利なのです。

スクリプトファイルをひとつ書いただけでの感想: 対話的・発見的な利用と、キチンとしたパッケージ作成のサポートは手厚いですが、少数のスクリプトファイルからなるプログラムの構造化・組織化が難しい印象です。

以下にソースファイル:

#
# ddds.R : Discrete-Discrete Dynamical System
# ###########################################

#
# 質点のぼかし処理
# =================

# (a, b)に置かれた質点を、パラメータvari(分散 = 2σ^2)でぼかした分布関数
blur <- function(vari, a, b, x, y) exp(-((x - a)^2 + (y - b)^2) / vari)

# (a, b)に置かれた質点をぼかした分布関数を返す
blurredPoint <- function(vari, a, b) function(x, y) blur(vari, a, b, x, y)

# 質点系に対するぼかした分布関数のz値マトリックスを作成する
makeBlurredPointsZMatrix <- function(
  # 質点系(複数の点のx座標, y座標)
  # ptXとptYは同じ長さのベクトル
  ptX, ptY,
  # それぞれの点の値(質量)
  vals,
  # ぼかしのパラメータ
  vari,
  # x-y平面のメッシュの定義
  meshX,  meshY
  )
{
  stopifnot(length(ptX) == length(ptY))
  zlist <- list()
  for (i in 1:length(ptX))
    zlist[[i]] <- vals[i]*outer(meshX, meshY, blurredPoint(vari, ptX[i], ptY[i]))
  Reduce(`+`, zlist, 0)
}

#
# dddsのオプション
# ================

# オプションを表示/設定する
ddds.options <- function(...) {
  if (!exists(".ddds.Options")) {
    .ddds.Options <<-
      list(
        # グラフのノード数
        nNodes = 10,
        # ぼかしのパラメータ
        vari = 0.01,
        # メッシュ生成のとき、区間を分割する小区間の数
        nSegments = 25,
        # 近隣の点として認識する距離
        neighborMax = 0.3,
        # 近すぎる点を除外するときの距離
        neighborMin = 0
        )
  }

  arglist <- list(...)
  opts <- .ddds.Options
  if (!length(arglist)) return(opts)

  if (is.null(names(arglist)) && !is.list(arglist[[1]])) {
    # 引数名はなく、第一引数の値があるとき
    # 引数は、問い合わせするオプション名の並び
    arglist <- unlist(arglist)
    if (length(arglist) == 1) opts[[arglist]]
    else opts[arglist]
  } else {
    # 引数名があるか、または第一引数がリスト
    origOpts <- opts
    if (is.list(arglist[[1]])) arglist <- arglist[[1]]
    # 設定を実行
    if (length(arglist) > 0) {
      opts[names(arglist)] <- arglist
      .ddds.Options <<- opts
    }
    invisible(origOpts)
  }
}

# オプションをリセット(再初期化)する
ddds.resetOptions <- function() {
  if (exists(".ddds.Options")) rm(".ddds.Options", inherits = TRUE)
  ddds.options()
}

#
# dddsオブジェクト
# ==================

# コンストラクタ
ddds <- function (
  # グラフの頂点の個数
  nNodes = ddds.options("nNodes"),
  # 頂点の位置
  x = NULL, y = NULL,
  # 隣接行列
  adj = NULL,
  # 遷移行列
  trans = NULL,
  # 以下、オプション値
  vari = ddds.options("vari"),
  nSegments = ddds.options("nSegments"),
  neighborMax = ddds.options("neighborMax"),
  neighborMin = ddds.options("neighborMin")
  )
{
  obj <- list()
  obj$nNodes <- nNodes
  obj$x      <- if (is.null(x)) runif(nNodes) else x
  obj$y      <- if (is.null(x)) runif(nNodes) else y
  obj$vari   <-  vari
  obj$nSegments   <- nSegments
  obj$neighborMax <- neighborMax
  obj$neighborMin <- neighborMin

  obj$vals   <- .ddds.genVals(obj$nNodes, obj$x, obj$y)
  obj$adj    <- 
    if (is.null(adj)) {
      .ddds.genAdj(obj$nNodes, obj$x, obj$y, obj$neighborMax, obj$neighborMin)
    } else {
      adj
    }
  obj$trans  <- if (is.null(trans)) .ddds.genTrans(obj$nNodes, obj$adj) else tarans
  
  structure(obj, class = "ddds")
}

#
# 内部使用関数
# ============

# 隣接行列を生成する
.ddds.genAdj <- function(
  nNodes, x, y,
  neighborMax = ddds.options("neighborMax"),
  neighborMin = ddds.options("neighborMin")
  )
{
  stopifnot(length(x) == nNodes, length(y) == nNodes)
  matVals <- rep(NA, nNodes^2)
  for (i in 1:nNodes) {
    matVals[(nNodes * (i -1) + 1) : (nNodes * i)] <- 
      as.numeric(
        y - y[i] >= 0
        & neighborMin^2 <= (x - x[i])^2 + (y - y[i])^2
        & (x - x[i])^2 + (y - y[i])^2 <= neighborMax^2
        )
  }
  matrix(matVals, c(nNodes, nNodes))
}

# 遷移行列を生成する
.ddds.genTrans <- function (nNodes, adj) {
  weights <- rep(NA, nNodes)
  for (j in 1:nNodes)
    weights[j] <- sum(adj[,j])
  trans <- matrix(NA, nNodes, nNodes)
  for (j in 1:nNodes) {
    if (weights[j] == 0) {
      trans[, j]  <- 0
      trans[j, j] <- 1
    } else {
      trans[, j] <- adj[, j]/weights[j]
    }
  }
  trans
}

# 頂点の初期値(質量分布)を生成する
.ddds.genVals <- function(nNodes, x, y) {
  vals <- rep(0, nNodes)
  k <- (which(min(y) == y))[1]
  vals[k] <- 1
  vals
}

# プロット用のメッシュの基準値を返す
mesh.ddds <- function(obj) {
  seq(0, 1, length = obj$nSegments + 1)
}

# プロット用にz値の行列を作成する
makeZMatrix.ddds <- function(obj) {
  ptX <- obj$x
  ptY <- obj$y
  vals <- obj$vals
  vari <- obj$vari
  meshX <- mesh.ddds(obj)
  meshY <- meshX
  makeBlurredPointsZMatrix(ptX, ptY, vals, vari, meshX, meshY)
}

#
# 外部へのインターフェース関数
# ============================

# 再初期化
Reinit.ddds <- function(
  obj,
  # 以下、オプション値
  vari = NULL,
  nSegments = NULL,
  neighborMax = NULL,
  neighborMin = NULL
  )
{
  obj$vari   <- if (!is.null(vari))     vari else obj$vari
  obj$nSegments   <- if (!is.null(nSegments))     nSegments else obj$nSegments
  obj$neighborMax <- if (!is.null(neighborMax)) neighborMax else obj$neighborMax
  obj$neighborMin <- if (!is.null(neighborMin)) neighborMin else obj$neighborMin
  
  obj$vals   <- .ddds.genVals(obj$nNodes, obj$x, obj$y)
  obj$adj    <- .ddds.genAdj(obj$nNodes, obj$x, obj$y, obj$neighborMax, obj$neighborMin)
  obj$trans  <- .ddds.genTrans(obj$nNodes, obj$adj)
  
  obj
}

# 状態遷移を行う
DoTrans.ddds <- function(obj, n = 1) {
  vals  <- obj$vals
  trans <- obj$trans
  for (i in 1:n) {
    vals <- trans %*% vals
  }
  obj$vals <- as.vector(vals)
  obj
}

# グラフィックデバイスに対してアニメーションを表示する
animate.ddds <- function (obj, n = 10, drawFun = persp, ...)  {
  stopifnot(n <= ani.options("nmax"))
  objlist <- list(obj)
  if (n >= 2) {
    for (i in 2:n) {
      obj <- DoTrans.ddds(obj)
      objlist[[i]] <- obj
    }
  }
  for (i in 1:n) {
    dev.hold()
    drawFun(objlist[[i]], ...)
    ani.pause()
  }
}

# オブジェクトをファイルにセーブする
save.ddds <- function(obj, filename = NULL) {
  filename <-
    if (!is.null(filename)) {
      filename
    } else {
      format(Sys.time(), "dddsObject-%Y%m%d_%H%M%S.RData")
    }
  message("saving object to '", filename, "'")
  save(obj, file = filename)
}

# ビデオアニメーションを生成する
genVideoAnimation.ddds <- function(obj, n = 10, drawFun = persp, ...) {
  saveVideo(
    {
      # アニメーション実行コード
      # ------------------------
      ani.options(interval = 0.1, nmax = 30)
      par(mar = c(3, 3, 2, 0.5),
          mgp = c(2, 0.5, 0),
          tcl = -0.3,
          cex.axis = 0.8, cex.lab = 0.8, cex.main = 1)
      animate.ddds(obj, n, drawFun, ...)
    },
    video.name = "ddds-example.mp4",
    # 以下、絶対パスが埋め込まれているのは好ましくない
    ffmpeg = "C:/Installed/ffmpeg/bin/ffmpeg.exe"
    )
}

#
# 汎用関数のメソッド
# ==================

plot.ddds <- function(obj, arrow = FALSE, ...) {
  n <- obj$nNodes
  x <- obj$x
  y <- obj$y
  adj <- obj$adj
  fun <- if(arrow) arrows else segments
  plot.default(x, y, xlim = c(0, 1), ylim = c(0, 1)) # グラフィックスの初期化
  for (i in 1:n)
    for (j in 1:n)
      if (adj[i,j]) fun(x[j], y[j], x[i], y[i], ...)
}

image.ddds <- function(obj, ...) {
  meshX <- mesh.ddds(obj)
  meshY <- meshX
  z <- makeZMatrix.ddds(obj)
  image.default(meshX, meshY, z, 
                xlim = c(0, 1), ylim = c(0, 1),
                xlab = "x", ylab = "y",
                ...)
}

persp.ddds <- function(obj,
                       theta = NULL,
                       phi = NULL,
                       expand = NULL,
                       col = NULL,
                       ...)
{
  theta  <- if (!is.null(theta)) theta   else 30
  phi    <- if (!is.null(phi)) phi       else 30
  expand <- if (!is.null(expand)) expand else 0.5
  col    <- if (!is.null(col)) col       else "lightblue"
  
  meshX <- mesh.ddds(obj)
  meshY <- meshX
  z <- makeZMatrix.ddds(obj)
  persp(meshX, meshY, z, 
        xlab = "x", ylab = "y", zlab = "z",
        xlim = c(0, 1), ylim = c(0, 1), zlim = c(0, 1),
        theta = theta, phi = phi, expand = expand, col = col,
        ...)
}

#
# 名前の管理
# ==========

# この“モジュール”で使用されている名前を列挙する
ddds.names <- function() {
  list(
    global = c("blur", "blurredPoint", "makeBlurredPointsZMatrix"),
    module = ls(envir = .GlobalEnv, pattern = "^ddds\\.*"),
    modulePrivate = ls(envir = .GlobalEnv, pattern = "^\\.ddds\\.*", all = TRUE),
    class = ls(envir = .GlobalEnv, pattern = "^.+\\.ddds$"),
    classPrivate = ls(envir = .GlobalEnv, pattern = "^\\..+\\.ddds$", all = TRUE)
    )
}

#
# 初期化の実行
# ============

if (!require('animation')) {
   install.packages('animation')
   stopifnot(require('animation'))
}
if (!exists(".ddds.Options")) ddds.resetOptions()