巨大なテーブルから SELECT する際のメモ

python経由で MySQL, PostgreSQL のメモリが爆発するレベルの巨大なテーブルからデータを取得する際のメモ(一千万件以上のレコードがあるテーブルから select するようなケースを想定)

現象

fetchall() を使わずに fetchone() / fetchmany() を使っているのに out of memory 的なエラーでプロセスが落ちる。
この現象は python の MySQLdb ライブラリ、および psycopg2 ライブラリで確認した。なお MySQL の場合、import MySQLdb ではなく、公式の import mysql.connector の場合はこの現象は生じない(のでこちらを使うべき?)。

原因

デフォルトのカーソルがクライアントサイドカーソルになっている。なので cursor.execute(sql) した時点ですべてのデータをクライアントに持ってくるので out of memory が発生する(このケースでは fetchone, fetchmany は「一つ(複数)ずつデータを取得する」という挙動をエミュレートしているに過ぎない)。

対策

サーバサイドカーソルを使う。やり方は以下のとおり。

MySQLdb

まずは失敗するパターン。

import MySQLdb
connection=MySQLdb.connect(
    host="host",user="user",
    passwd="passwd",db="db")
cursor=connection.cursor()
cursor.execute("select * from very_huge_table")  # 失敗!
row = cursor.fetchone()
while row is not None:
    print row
    row = cursor.fetchone()

SSCursor を使えばうまくいく。

import MySQLdb
import MySQLdb.cursors   # 追加!
connection=MySQLdb.connect(
    host="host",user="user",
    passwd="passwd",db="db",
    cursorclass = MySQLdb.cursors.SSCursor) # 追加!
cursor=connection.cursor()
cursor.execute("select * from very_huge_table")
row = cursor.fetchone()
while row is not None:
    print row
    row = cursor.fetchone()

psycopg2

psycopg2 の場合はカーソルオブジェクトを生成する際に name を指定するとサーバサイドカーソルになる、という記述があった。

参考URL : http://initd.org/psycopg/docs/connection.html#connection.cursor

If name is specified, the returned cursor will be a server side cursor (also known as named cursor). Otherwise it will be a regular client side cursor.

以下のコードは execute の時点でクライアントサイドにすべてのデータを引っ張ってこようとして失敗する。

import psycopg2
connection = psycopg2.connect("dbname=db host=host user=user")
cursor=connection.cursor()
cursor.execute("select * from very_huge_table")  # 失敗!
row = cursor.fetchone()
while row is not None:
    print row
    row = cursor.fetchone()

named cursor を使えばOK。

import psycopg2
connection = psycopg2.connect("dbname=db host=host user=user")
cursor=connection.cursor(name="thecursor") # 違いはここだけ。名前は PostgreSQL の非予約語ならなんでも
cursor.execute("select * from very_huge_table")
row = cursor.fetchone()
while row is not None:
    print row
    row = cursor.fetchone()

image.plot の目盛を axis 関数で描く

Rのfieldsパッケージのimage.plot関数で目盛に数値ではなくテキストを書き込みたい。axis関数を使っても何故か反映されないというバグ?がある(fieldsのバージョンは6.9.1)。例えば以下の例:

library(fields)

x <- matrix(rnorm(25),nc=5)

colnames <- paste0("col",1:5)
rownames <- paste0("raw",1:5)

# Fail
image.plot(1:5, 1:5, x, xaxt="n", yaxt="n", xlab="X-axis", ylab="Y-axis")
axis(side=2,at=1:5,labels=rownames)
axis(side=1,at=1:5,labels=colnames)

image_plot_fail

確かに最後の二行の axis 関数の結果が反映されていない(x軸, y軸の目盛に文字が書けていない)。

r- how to edit elements on x axis in image.plot
この問題はここでも議論されていて、普通のimageを使って、凡例だけimage.plotを使え、とある。

もう少し簡単に解決する方法を見つけた。

# Success
image.plot(1:5, 1:5, x, xaxt="n", yaxt="n", xlab="X-axis", ylab="Y-axis")
text(-100,-100,".")  # This is dummy, but required for drawing axis
axis(side=2,at=1:5,labels=rownames)
axis(side=1,at=1:5,labels=colnames)

image.plot と axis の間にダミーでtext関数を実行する(三行目)。以下のとおり不思議なことにこれで解決した(理由は未調査)。

image_plot_success

Ruby CGI で Insecure operation error

久しぶりの ruby ネタ。以下のような ruby CGI で 500 internal error が出た。

#!/usr/bin/ruby

require "cgi"
cgi = CGI.new
params = cgi.params

open(params["filename"],"w") { |f| f.puts "hoge" }
cgi.out { "hoge" }

apache のログを見てみるとファイルを作成する部分 (open) で Insecure operation というエラーが出ている。

いろいろ調べてみると、どうやら ruby CGI で URL のパラメータを使ったパス名でファイルを書き込もうとしたことが原因のようだ。パス名に空白とか入っていると悪さができそう、ということらしい。ごもっともです。
“Ruby CGI で Insecure operation error”の続きを読む

[R] データフレームを可視化するtabplotパッケージについて

データの列数が多くなってくる (高次元になってくる) とデータの全体像が捉えにくくなる.R-bloggers を読んでたら,Rのデータフレームを可視化するためのパッケージ tabplot なんてのがあるらしい,ということで試してみた.コレを使うとデータを効率的に可視化できるので,よくわからないデータをとりあえず可視化してみて,それからあれこれ考えてみると捗るかもしれない.

CRAN から tabplot パッケージをインストールすればすぐに使える (類似の名前で tableplot なんてのがあるので注意).
“[R] データフレームを可視化するtabplotパッケージについて”の続きを読む

Google検索時に出てくるevernoteのアレの位置を変更する

Chrome で evernote プラグインを使ってるときに google 検索すると evernote の検索結果も同時に出てくるアレの話です。ピンとこない人は関係ないはずです。

あの検索結果は個人的には結構便利だとは思うのですが、レイテンシがわりと大きいため、一番上のリンクをクリックしようとした瞬間にあの水色のボックスが出現したりすることが多いのが少し困る。

かといって機能を切ってしまうのもなんとなく忍びないのですこし改良してみた。最初は javascript だけでやろうとしたけど、手持ちの chrome には stylebot という拡張機能がついていたのでそれを使うことにした。
“Google検索時に出てくるevernoteのアレの位置を変更する”の続きを読む

Rで階段プロットを書く

知らなかったのでメモ。Rで階段プロットを書くには type=”s” を指定する。

以下は経験累積分布関数を書くためのコード片。

cumplot <- function(x,...){
  plot(sort(x),(1:length(x))/length(x),type="s",...)
}

cumplot(rnorm(100),xlab="",ylab="",xlim=c(-3,3))
curve(pnorm,add=T,col=2)
legend("topleft",legend=c("empirical CDF","normal CDF"),
       lt=1,col=1:2)

こんな図が書ける。

ほかにも type=”S” (大文字) なんてのがあって、だいたい同じなんだけど、S は上にずれる。経験累積分布の場合は小文字 s が正解。

Rにおける値渡しと参照渡し (2)

昨日の記事を書いたあと、twitter で tracemem 使うといいよ、と教えていただきました (@sfchaosさん、ありがとうございました!)。

この関数ははじめて知ったのですが、ヘルプを見て意訳すると以下の様な感じです。

tracemem(x) を実行すると x について duplicate (複製)が生じた時にメッセージを表示する。これは2つのオブジェクトがメモリを共有している場合において、1つが変更された時に生じる。ちなみに untracemem(x) でメッセージフックを解除できる。

ということで関数の内部でコピーが起きたか、を判別するのにピッタリです。

今回の内容:

  • tracemem を使って値渡し、参照渡しを確かめる
  • 参照渡しなのか、値渡しなのか、が確定するタイミングについて考える

tracemem を使った方法

以下の方法は @sfchaos さんに教えてもらった方法そのままです。前回同様、prod1 は行列 A について値渡しになりそうな関数で、prod2 は参照渡しになりそうな関数です。

prod1 <- function(A,x){
  A <- A + diag(x)
  A %*% x
}

prod2 <- function(A,x){
  A %*% x
}

N <- 1000
A <- matrix(rnorm(N*N),nc=N)
x <- rnorm(N)

tracemem(A)
# => [1] "<0x7e790008>"
invisible(prod1(A,x))  # invisible は結果を print しないようにする関数
# => tracemem[0x7e790008 -> 0x7dfe0008]: prod1  # コピーが発生した!

invisible(prod2(A,x))
# => 何も表示されない!=コピーが生じていない!

ということで、前回は実行時間から推論しただけでしたが、やはり引数を変更するとコピーが生じる、ということで間違いないことが確認できました。

複製はどのタイミングで生じるのか?

前回、以下のように書きました。

そして、引数が変更されるかされないかはパースした段階でわかる(なのでパースの段階で値渡しか参照渡しかを判別することが可能)

ですが、これは誤りでした。実際、以下のような関数はパースの段階で値渡しか参照渡しかを判別できるでしょうか?

prod3 <- function(A,x,add=T){
  if(add) A <- A + diag(x)
  A %*% x
}

add が真のときは引数が変更されて、偽のときは引数が変更されません。

これをパースの時点で判別しようとすると、Rの副作用を作らないという原則から add がいかなる値であろうともコピーを行うことになるはずです。

もう一つの可能性としては、実際に変更が起きたその瞬間にコピーを作る、というものがありえるでしょう。

実験してみます。

N <- 1000
A <- matrix(rnorm(N*N),nc=N)
x <- rnorm(N)

tracemem(A)
# => [1] "<0x7dfe0008>"

invisible(prod3(A,x,add=T))
# => tracemem[0x7dfe0008 -> 0x7ef40008]: prod3

invisible(prod3(A,x,add=F))
# => 何も表示されない!=コピーが生じていない!

同じ関数でも引数の状態によってコピーが発生したり、しなかったり。add=F のときは prod3 の if の内部まで進まないため、A が変更されず、したがってコピーが発動しない、ということになります。

つまりRは

  • 引数のコピーを作るか(値渡しか)、作らないか(参照渡しか)は実際に引数が変更される瞬間ギリギリまで判別しない、
  • 変更される瞬間(直前?)で重い腰を上げてコピーを作成する(遅延評価)

という挙動をしているようですね。

まとめ

  • Rは基本的には値渡しであり、C++のように引数として与えられた変数を変更することで外側の世界に影響させることはできない
  • 関数内部で引数を変更しない場合は、コピーが生じない (C++ の const& のようなイメージ)
  • 関数内部に引数を変更するコードがあったとしても、実際に引数が変更される段まではコピーは生じない (遅延評価)

Rにおける値渡しと参照渡し

Rの関数に引数を渡すと値渡しになる、とずっと信じていたわけだけど、どうも違うらしい。Rはどうやら「自動的に」値渡しすべきか、参照渡しにすべきか、を判断しているようだ。

C++みたいな言語では「フツーに」引数を渡すと全部値渡しになって(特に行列のような巨大なオブジェクトを渡す場合には)効率が悪いので、ポインタや参照で渡したり、副作用を気にする場合は const 参照で渡したりする。

さて、上で述べた R の自動判断機能は以下のような事実に基づくようだ。

  • そもそも値渡しは関数に引数を渡した段階でオブジェクトのコピーを生成して「引数として渡された変数を関数の内部で変更してもスコープの外では値が変更されない」ことを保証するのだが、
  • オブジェクトが関数の内部で変更されないことが保証されるならば「関数の実行中はスコープの外でも値が変更されない」ことは保証されるのでコピーをそもそも生成する必要がない、
  • そして、引数が変更されるかされないかはパースした段階でわかる(なのでパースの段階で値渡しか参照渡しかを判別することが可能) → 間違いでした。詳細はRにおける値渡しと参照渡し(2)に書きました。

実験

これは以下のようにして確かめることができる、はず。

prod1 <- function(A,x){
  A[1,1] <- A[1,1] + 1
  A %*% x
}

prod2 <- function(A,x){
  x[1] <- x[1] + 1
  A %*% x
}
  • prod1, prod2 はともに引数として渡された行列 A とベクトル x の積を計算する
  • prod1 では A[1,1] に 1 を加える
  • prod2 では x[1] に 1 を加える(これはprod1の代入作業のコストと揃えるため)
  • その後、積を計算する

ということをやっているのだが、上で述べたようなことがただしければ、

  • prod1 では行列 A のコピーが発生し
  • prod2 ではベクトル x のコピーが発生する

ため、より巨大なオブジェクトを渡すことになる prod1 が (生じる四則演算の数は同一にもかかわらず) 大幅に遅くなることが予測される。

実際、1000 x 1000 程度の行列を考えると以下の様な結果が得られる。


N <- 1000
A <- matrix(rnorm(N*N),nc=N)
x <- rnorm(N)

system.time( for(i in 1:1000) prod1(A,x) )
# =>   user  system elapsed
# =>  11.03    3.50   14.56

system.time( for(i in 1:1000) prod2(A,x) )
# =>  user  system elapsed
# =>  4.59    0.01    4.62

prod1のほうがかなり遅くなっていることが分かる。これより、上で述べたことはおそらく正しいだろうと結論される。

まとめ

このことを知ったからといって高速なプログラムが書けるようになるわけではないが、少なくとも「こんな大きなオブジェクトを関数に渡すのは気が引ける(なのでグローバル変数にしてしまえ)」という不安を払拭することができると思う。

【追記】続きを書きました。

Rからc++コードを呼び出すRcpp+inlineを試す

従来、RをCで拡張するときは「Rの基礎とプログラミング技法(→amazon)」にもあるとおり、R.h をインクルードして SEXP (S EXPression) 型と格闘するイメージだけど(何度か実践で使ったけど気楽にやろう、ってものではなかった)、Rcpp + inline を使うと魔法のように簡単に R と c++ を連携させることができるようだ。

なんとインラインで c++ を書いて、それを直接実行する、というすごい実装(inlineパッケージを使った場合。Rcpp は R.h のラッパー)。

使い道としてはMCMCとかやたら遅いアルゴリズムを実装するときに便利そう(C++の中でRの豊富な乱数生成が使えたりするから)。

他にもループが遅い時にループの深いところを c++ で書き換えてしまうとか。かなり気楽だからこそできるワザ。

勉強を兼ねて適当なテスト関数を書いてみた。
“Rからc++コードを呼び出すRcpp+inlineを試す”の続きを読む

Rの配列/リストのランダムアクセスのパフォーマンスを調べてみた

よくRでは配列/リストに対してfor文+インデックスでアクセスしてはいけない、ということが言われますが、そうせざるを得ない場面にも時々遭遇します。少し気になってインデックスでアクセスした時のパフォーマンスの特性を調べて見ました。

通常、データ構造の意味で配列といえばランダムアクセスは定数時間、つまり O(1) で、リストのランダムアクセスは O(n) です。ところが、Rの場合は必ずしもそうなっていないようです

実験条件

2つのアクセスパターンで実験をしました。

■ アクセスパターン1
配列とリストの i 番目の成分に100万回 1 を加える。i=1,2000,4000,6000,…,20000 に対して実行時間を調べる。

■ アクセスパターン2
配列とリストの 1~i 番目の中のランダムに選んだ添字の成分に 1 を加える操作を100万回行う。i=1,2000,4000,6000,…,20000 に対して実行時間を調べる。

結果

2つのアクセスパターンに対する結果を表したのが次の図(使ったスクリプトは記事の一番最後):

縦軸は100万回のランダムアクセスに要した時間、横軸は上の実験条件の所で述べた i を示しています。

考察

はっきり言ってうまく実験がうまく制御されているかどうかについては自信がありませんが(定数分のオーバーヘッドがどこに起因しているのかよくわからない)、この実験の結果を真とすればアクセスパターン1では配列/リストともに定数時間であり、一方、アクセスパターン2では配列は定数時間ですが、リストの方は長さが増加するとアクセス時間が増加することになります。
この実験をする前は

  • リストと名のつくものは無条件で所要時間が線形に増えるだろう

と思っていたので少し新鮮な結果でした。

疑問・まとめ

  • リストの同じ添字にアクセスするのは定数時間なのは、最近アクセスした要素を何らかの方法でキャッシュしているから?
  • 配列に関して、パターン1よりもパターン2のほうが所要時間が短い理由はなぜ?

まあ、100万回程度のアクセスにこんなに時間がかかるならやっぱりできるだけ使わないほうがよいということでしょう(100万次元のベクトルの足し算などはこの実験よりもはるかに高速に行うことができる)。どうしても使う場合は、特にリストの場合はキャッシュが効くように同じ添字になるべく連続してアクセスすべし、ということを頭の片隅にいれておくことにします。

実験に使ったスクリプト

ntry <- 1000000 # Number of repeat

ary.len <- 20000

# Defining index for test
idx <- seq(0,ary.len,ary.len/10)
idx[1] <- 1

# Array access time --- Pattern 1
ary <- rep(0,ary.len)
durations.ary <- matrix(0,length(idx),5)
j <- 1
for(i in idx){
  n <- sample(1:i,ntry,T)
  start <- proc.time()
  for(k in 1:ntry)
    ary[i] <- ary[i] + 1
  durations.ary[j,] <- as.numeric(proc.time() - start)
  j <- j+1
  gc();gc();gc();
}

# List access time --- Pattern 1
lst <- list(0)
for(i in 1:ary.len)
  lst[[i]] <- 0
durations.lst <- matrix(0,length(idx),5)
j <- 1
for(i in idx){
  n <- sample(1:i,ntry,T)
  start <- proc.time()
  for(k in 1:ntry)
    lst[[i]] <- lst[[i]] + 1
  durations.lst[j,] <- as.numeric(proc.time() - start)
  j <- j+1
  gc();gc();gc();
}

tmp <- c(durations.ary[,3],durations.lst[,3])
par(mfrow=c(1,2))
plot(idx,durations.ary[,3],type="n",ylim=c(7,8),xlab="Index of array/list",ylab="Excution time")
lines(idx,durations.ary[,3],type="o",lw=2,col="steelblue")
lines(idx,durations.lst[,3],type="o",lw=2,col="pink")
legend("bottomright",legend=c("Array access time","List access time"),col=c("steelblue","pink"),pch=1,lw=2)
title("Access pattern 1")


# Array access time --- Pattern 2
ary <- rep(0,ary.len)
durations.ary2 <- matrix(0,length(idx),5)
j <- 1
for(i in idx){
  n <- sample(1:i,ntry,T)
  start <- proc.time()
  for(m in n){
    ary[m] <- ary[m] + 1
  }
  durations.ary2[j,] <- as.numeric(proc.time() - start)
  j <- j+1
  gc();gc();gc();
}

# List access time --- Pattern 2
lst <- list(0)
for(i in 1:ary.len)
  lst[[i]] <- 0
durations.lst2 <- matrix(0,length(idx),5)
j <- 1
for(i in idx){
  n <- sample(1:i,ntry,T)
  start <- proc.time()
  for(m in n){
    lst[[m]] <- lst[[m]] + 1
  }
  durations.lst2[j,] <- as.numeric(proc.time() - start)
  j <- j+1
  gc();gc();gc();
}


tmp <- c(durations.ary2[,3],durations.lst2[,3])
plot(idx,durations.ary2[,3],type="n",ylim=c(7,8),xlab="Index of array/list",ylab="Excution time")
lines(idx,durations.ary2[,3],type="o",lw=2,col="steelblue")
lines(idx,durations.lst2[,3],type="o",lw=2,col="pink")
legend("bottomright",legend=c("Array access time","List access time"),col=c("steelblue","pink"),pch=1,lw=2)
title("Access pattern 2")