RからC言語の関数を呼ぼう
Rから呼ばれるためのC言語の関数は、高速で並列処理がOpenMPで容易でRでサポートされないOSのAPIやライブラリを使えるだけではなく、Rの機能も大きな制限なく操作することができます。
- 1 準備
- 2
.Cと.Callと.Externalの3種類のC呼び出し - 3
.Cを使ってみよう - 4
.Callを使ってみよう - 5 DLLで引数の型や呼び出し名を指定する
- 6
可変長引数のときに使う
.External - 7 C言語から利用できるRのAPI
- 8 OpenMPによる並列化
- 9 Rのゾンビがメモリに残ってファイルを放さない！
- 10 まとめ
統計解析用プログラミング言語のRの長所の一つは、C言語で機能を拡張しやすいことです。単にCで書かれた関数を呼べるだけではなく1、CにRのオブジェクトを操作する機能と、数学確率関係のビルトイン関数を提供します。また、ボイラープレートコード2が少ないです。
特殊なライブラリをリンクするCコードのコンパイルオプションは煩雑になり勝ちですが、Cコードを用意する簡単なコマンドが用意されており、その苦労もありません。WindowsではCコンパイラなどを用意するのが苦労なときもありますが、Rtoolsと言う開発環境が提供されています。
CによるR拡張の応用にRcppがあり、主に計算速度の向上を目指す場合はRcppを使う方が少ない労力3で可読性が高いコードを生産できますが、その土台になっている部分だけでも有用です。
ドキュメントも何だかんだとある方で、最新のWriting R Extensions4をよく読み、インクルードファイルを追いかけたら使い方で分からないことは無いと思いますが、最初はとっつきづらい面もあるので、簡単な例を紹介していきたいと思います。
1 準備
MS-WindowsではRとRtoolsをインストールして、環境変数が適切に設定されているか確認してください。この文書の記述時点での最新のR
4.2.2とrtools
4.2の標準インストールであれば、環境変数RTOOLS42_HOMEがC:\rtools42になります。また、コマンドプロンプトからRを呼び出せるように、PATHに%R_HOME%\binを追加しておいてください5。
LinuxであればOSのパッケージマネージャーを用いるなどしてgfortranとRが動く環境にしてください。
インストールが済んだら、OpenMPが有効になるようにオプションを設定します。Windowsの場合は%R_USER%\.R\Makevars.winを作成します。Linuxの場合は~/.R/Makevarsです。Rcppなどで開発がしたことがある人にはお馴染みだと思います。
Windowsでこのファイルが存在しない場合は、Rで
directory <- file.path(Sys.getenv("R_USER"), ".R")
if(!dir.exists(directory)) dir.create(directory)
cat("CFLAGS = $(SHLIB_OPENMP_CFLAGS)\nFCFLAGS = $(SHLIB_OPENMP_FFLAGS)",
file=paste(directory, "Makevars.win", sep="/"))とすれば作られます。場所はprint(directory)をして確認してください。
作ったMakevars.winの中身のCFLAGSはCコンパイラに対するオプション、FCFLAGSはgFortran（のfree
source
formのコンパイル）に対するオプションになります。それぞれ$(SHLIB_OPENMP_CFLAGS)と$(SHLIB_OPENMP_FFLAGS)を指定しておくと、OpenMPが有効になります。Cコンパイラは不要に思えるかもですが、gFortranがつくったオブジェクト・ファイルをリンクするのに使います。
2
.Cと.Callと.Externalの3種類のC呼び出し
Rには.Cと.Callと.Externalの3種類のC呼び出しがあります。.CはただのCの関数を呼ぶI/Fで、.Callと.ExternalはRのために書かれたCの関数を呼ぶI/Fです。Cの関数の引数や戻り値をRのオブジェクトにする場合は、.Callか.Externalを使うことになります。.Callは固定長呼び出し、.Externalは可変長呼び出しと言う違いがあります。.Cは古い規格で、現在では.Callか.Externalを使うことが推奨されています。実際、Rらしい挙動の関数を書く場合は、.Callか.Externalを使うしかないです。
3
.Cを使ってみよう
以下の内容のexample_add.cと言うファイルを用意し、
以下のようにコンパイルを行い、
R CMD SHLIB example_add.c Rの作業ディレクトリをexample_add.cにあるフォルダーに移動すれば、以下のようにRから呼び出すことができます。
dll <- paste("example_add", .Platform$dynlib.ext, sep="")
dyn.load(dll)
.C("add", as.double(123), as.double(321))
dyn.unload(dll)[[1]]
[1] 444
[[2]]
[1] 321リストの最初の要素が引数123と321の合計の444になっているので、Cのコードが動いたことが分かります。
DLLのファイルパス名をつくってDynamic Link Library（以下、DLL）をロードし、C関数を呼んだ後に（もうC関数が不要になったので）DLLをアンロードしているだけです。引数の型はCの関数にあわせて変換しておきます。
Cの関数はvoid型で戻り値はなく、引数のポインターの先を更新することで、計算結果を戻します。ベクトルを引数として渡せますが、ベクトルの長さは別の引数で指定する必要があります。具体的に使えるRの引数の型とC言語における対応は、Writing
R Extensionsの以下の表にまとまっています。
| R storage mode | C type | Fortran type |
|---|---|---|
| logical | int * | INTEGER |
| integer | int * | INTEGER |
| double | double * | DOUBLE PRECISION |
| complex | Rcomplex * | DOUBLE COMPLEX |
| character | char ** | CHARACTER(255) |
| raw | unsigned char * | none |
Rcomplexを使うときは、Cのソースコードに#include <R.h>が要ります。characterはポインタのポインタのcharになるので、引数と異なるサイズの文字列に置き換えることができますが、他は引数のサイズに制約されます。
4
.Callを使ってみよう
.Cの機能で十分なことも多いと思いますが、引数と戻り値の型が制約されているので、現在では.Callと.Externalの利用が推奨されています。
.Callと.Externalは引数と戻り値がSEXP型になっています。Rのオブジェクトは内部的にSEXPREC構造体に保存されているのですが、それへのポインターになる型です。つまり、listやvectorを引数にとり、戻り値に使うことが可能になります。
4.1 ベクターを引数にとって、ベクターを戻す関数
習うより慣れろ感があるので、実際にコードを見て行きましょう。まずは、ベクターを引数にとって、ベクターを戻す関数を作ってみます。
以下の内容のファイルexample_vector.cを作成し、
#include <R.h>
#include <Rinternals.h>
/* 引数も戻り値もSEXP型 */
SEXP abs_roundup(SEXP v)
{
int i;
double x;
SEXP r, names;
/* 型チェック */
if(!isVector(v))
error("'v' should be a vector.");
if(!isReal(v))
error("'v' should be of numeric.");
/*
戻り値になるベクターを作成
REALSXPは実数
length(v)でvの要素数を取得
allocVector(...)で、ベクトルを作成
PROTECT(...)で、ガーベッジコレクターから保護
*/
PROTECT(r = allocVector(REALSXP, length(v)));
for(i=0;i<length(v);i++){
/* 欠損値か判定（説明用で無くても動作は同じ） */
if(ISNA(REAL(v)[i])){
/*
REAL(r)は、double*のポインターを意味するマクロ
ポインタに慣れていないRユーザーからポインタを隠蔽する
*/
REAL(r)[i] = NA_REAL;
continue;
}
/* 例えば2.3は3に、-2.3は-3にする */
if(0>(x = REAL(v)[i]))
REAL(r)[i] = -1*ceill(fabs(x));
else
REAL(r)[i] = ceill(x);
}
/*
ベクターrの'names'属性として、ベクターvの'names'属性を代入する
なお、`names`属性のデータもSEXP型のベクター
*/
setAttrib(r, R_NamesSymbol, getAttrib(v, R_NamesSymbol));
/* GCからの保護を解除 */
UNPROTECT(1);
return r;
}以下の内容のexample_vector.Rを作成し、
# DLLのパス名を作ってロードする
dll <- paste("example_vector", .Platform$dynlib.ext, sep="")
dyn.load(dll)
# C関数に引数として与えるベクトルを作成
v <- c(-1, sqrt(2), -1*sqrt(3), sqrt(4), -1*sqrt(5), NA)
names(v) <- c("1st", "2nd", "3rd", "4th", "5th", "6th")
# 引数を確認
print("input:")
print(v)
# C関数を呼び出して結果を得る
print("output:")
print(.Call("abs_roundup", v))
# 比較のためceilingの結果を示す
print("ceiling:")
print(ceiling(v))
# DLLをアンロードする
dyn.unload(dll)2つのファイルのあるディレクトリで、
R CMD SHLIB example_vector.c
R --slave -f example_vector.Rとコンパイルと実行を行なうと、
[1] "input:"
1st 2nd 3rd 4th 5th 6th
-1.000000 1.414214 -1.732051 2.000000 -2.236068 NA
[1] "output:"
1st 2nd 3rd 4th 5th 6th
-1 2 -2 2 -3 NA
[1] "ceiling:"
1st 2nd 3rd 4th 5th 6th
-1 2 -1 2 -2 NAと結果を得ることができます。符号つきで切り上げるceilingと似て非なる、絶対値を切り上げて符号を戻す関数ができました。3rdと5thが異なりますね。
なお、sign(v)*ceiling(abs(v))で同じことができてかつ高速なので、このC関数自体は役立ちません。あくまでベクトルを扱う例なので悪しからず。
4.1.1
error関数
void error(const char * format, ...)はエラーを生じさせる関数で、void warning(const char * format, ...)は警告を生じさせる関数です。単に引数を出力するだけかと思いきや、何とprintfと同じ構文が取れる優れものです。ライブラリでは計算の前提があわないときはエラーにするしかないので有用ですし、printfと同等の機能があってエラーや警告の理由を詳しくメッセージに載せられるので重宝します。
4.1.2
SEXP型の種類（保持しているデータの型）
Rのベクターが数値や文字列といった型を持つように、RのSEXP型（の参照先のSEXPREC構造体）も内部に型を持ちます。上の例ではREALSXPのベクターを作っていました。
| Rの型 | SEXP型の種類 |
|---|---|
| numeric | REALSXP |
| integer | INTSXP |
| logical | LGLSXP |
| single | SINGLESXP |
| character | STRSXP |
| list | VECSXP |
Rのベクターと異なり、SEXP型はlistも表すことに注意してください。
型変換はnew_vector = PROTECT(coerceVector(old_vector, REALSXP));と言う風にできます6が、bad
practiceになり勝ちなので避けましょう。
4.1.3
is*関数
変数の型はTYPEOF(SEXP)を使ってINTSXPやREALSXPなどになっているか確認してもよいのですが、is*関数も用意されています。
Rboolean isVector(SEXP);
Rboolean isVectorAtomic(SEXP);
Rboolean isVectorList(SEXP);
Rboolean isMatrix(SEXP);
Rboolean isPairList(SEXP);
Rboolean isPrimitive(SEXP);
Rboolean isTs(SEXP);
Rboolean isNumeric(SEXP);
Rboolean isArray(SEXP);
Rboolean isFactor(SEXP);
Rboolean isObject(SEXP);
Rboolean isFunction(SEXP);
Rboolean isLanguage(SEXP);
Rboolean isNewList(SEXP);
Rboolean isList(SEXP);
Rboolean isOrdered(SEXP);
Rboolean isUnordered(SEXP);4.1.4 欠損値，非数値，無限
欠損値や無限のテストには以下のマクロが使えます。
| 関数 | R’s NA | IEEE NaN | Inf/-Inf |
|---|---|---|---|
| ISNA(x) | ✓ | ||
| ISNAN(x) | ✓ | ✓ | |
| R_FINITE(x) | ✓ | ✓ | ✓ |
実数値としてはNA_REAL，R_NaN,R_PosInf,R_NegInfが用意されています。
4.1.5 PROTECTとUNPROTECT
Rの内部の挙動をコントロールするおまじないのようなものですが、allocVectorしてベクトルを作ったあと、ベクトルの要素をセットするまではPROTECTマクロでガーベッジコレクターからSEXP型が指し示すSEXPREC構造体を保護します。スタック構造になっており、保護を解除するときはUNPROTECTマクロで解除する数を指定します。複雑なリスト構造の戻り値をつくったときなどは、UNPROTECT(4)のように複数同時に保護解除もできます。
4.1.6 属性の取得と付与
getAttrib関数で引数の変数からnames属性のベクトルを取得し、setAttribで戻り値になる変数にそれを設定しているわけですが、Rでの属性名とCのSymbol
Tableには以下のような対応があります。
| R | C |
|---|---|
class |
R_ClassSymbol |
names |
R_NamesSymbol |
dim |
R_DimSymbol |
dimnames |
R_DimNamesSymbol |
levels |
R_LevelsSymbol |
直接使う機会が多いのは以上だと思いますが、Rinternals.hを参照すると他のSymbolも分かります。
4.2 行列を引数にとって、行列を戻す関数
Rの行列はベクトルに属性がついただけの構造なので、ベクトルと同様に処理できます。ライフゲームの1回だけ動かす関数を作ってみましょう。ライフゲームは生命の誕生や死を計算機上でシミュレーションするゲームで、今回は典型的な周囲のセルに応じて生まれたり(=1)死んだり(=0)するものを作成します。
以下の内容のファイルLifeGame.cを作ります。
#include <R.h>
#include <Rinternals.h>
SEXP LifeGame(SEXP m)
{
SEXP ans, dim;
int nor, noc, i, j, nol, x, y, r, c;
/* 第一引数が行列か確認 */
if(!isMatrix(m)){
error("A matrix is required for the first argument.");
}
/* 型チェック*/
if(!isInteger(m)){
error("A integer matrix is required.");
}
/* 行列の行数、列数を得る */
dim = getAttrib(m, R_DimSymbol);
nor = INTEGER(dim)[0]; /* 行数 */
noc = INTEGER(dim)[1]; /* 列数 */
/* 戻り値の行列を作成（注意：REALSXPではなくINTSXPを指定） */
PROTECT(ans = allocMatrix(INTSXP, nor, noc));
/* 端は消去 */
for(i=0; i<nor; i++){
INTEGER(ans)[i] = 0;
INTEGER(ans)[i + nor*(noc-1)] = 0;
}
for(j=0; j<noc; j++){
INTEGER(ans)[nor*j] = 0;
INTEGER(ans)[nor-1 + nor*j] = 0;
}
/* 中央部分だけを処理 */
for(i=1; i<nor-1; i++){
for(j=1; j<noc-1; j++){
/* 周辺の8マスの状態を確認 */
nol = 0;
for(y=-1; y<=1; y++){
for(x=-1; x<=1; x++){
if(!y && !x)
continue;
if(0<INTEGER(m)[i+y + nor*(j+x)])
nol++;
}
}
/* 周囲のマスの状態で生死を決定 */
if(0<INTEGER(m)[i + nor*j])
INTEGER(ans)[i + nor*j] = nol<=1 ? 0 : (nol>=4 ? 0 : 1);
/* 周囲のマスの状態で発生を決定 */
else
INTEGER(ans)[i + nor*j] = nol==3 ? 1 : 0;
}
}
/* 引数のrownamesとcolnamesを戻り値にコピーする */
setAttrib(ans, R_DimNamesSymbol, getAttrib(m, R_DimNamesSymbol));
UNPROTECT(1);
return(ans);
}コンパイルします。
R CMD SHLIB LifeGame.c次に、LifeGame.cがあるディレクトリに、以下の内容のファイルLifeGame.Rを作成します。
# ライブラリを呼び出す
dll <-paste("LifeGame", .Platform$dynlib.ext, sep = "")
dyn.load(dll)
# Rの関数でラッピングして使う
lgame <- function(m){
.Call("LifeGame", m)
}
# グライダー
mg <-matrix(as.integer(c(
1,1,1,
1,0,0,
0,1,0)), 3, 3, byrow=TRUE)
# 移動物体なので広めの行列を作る
m <- matrix(as.integer(rep(0,12*12)),12,12)
# 広い行列にグライダーを写す
for(r in 1:nrow(mg)){
for(c in 1:ncol(mg)){
m[r+nrow(m)-nrow(mg)-1,c+ncol(m)-ncol(mg)-1] <- mg[r,c]
}
}
# 移動させてみる
for(i in 1:15){
print(m)
m <- lgame(m)
}
# DLLをアンロードする
dyn.unload(dll)二つのファイルがあるディレクトリで、
R --slave -f example_vector.Rと入力すると、結果をここに書くのは省略しますが、1と0で表現されたライフゲームの15ステップが連続で表示されます。|more（Linuxだと|less）をつけないと、ちょっと状況が把握できないかもですが。
4.3 リストを引数にとって、リストを戻す関数
Rの構造データはリストになるので、Cでリストを処理できるとデータの受け渡しがスムーズになります。
引数のリストの中の整数、実数、複素数のベクトルは合計値を、それ以外は引数のリストの中をコピーして戻す関数を作ってみましょう。
以下の内容のファイルexample_list.cを作ります。
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Complex.h>
/* 引数も戻り値もSEXP型 */
SEXP listsum(SEXP lst)
{
SEXP r, e, k;
int i, j, s_i;
double s_d;
Rcomplex in;
/* 型チェック/Listは内部的にはNewList */
if(!isNewList(lst))
error("A list is required for the first argument.");
/* 戻り値となるlength(lst)の大きさのリスト（になるベクター）を確保 */
PROTECT(r = allocVector(VECSXP, length(lst)));
for(i=0; i<length(lst); i++){
/* Listのi番目の要素を取り出す */
e = VECTOR_ELT(lst, i);
/* Integer/Real/Complexは入力ベクトルの合計値を、それ以外は入力ベクトルをコピー */
switch(TYPEOF(e)){
case INTSXP:
for(j=0,s_i=0; j<length(e); j++)
s_i += INTEGER(e)[j];
/* ScalarIntegerを使わないと長くなる */
/* 要素の数が1のベクターを作成, GCから保護 */
k = PROTECT(allocVector(INTSXP, 1));
/* 要素に値を代入 */
INTEGER(k)[0] = s_i;
/* Listのi番目の要素としてベクターをセット */
SET_VECTOR_ELT(r, i, k);
/* GCからの保護を解除 */
UNPROTECT(1);
break;
case REALSXP:
for(j=0,s_d=0; j<length(e); j++)
s_d += REAL(e)[j];
/* double型の引数を長さ1の配列の要素の値にしたSEXP構造体に変換するScalarRealを使うと短く済む */
SET_VECTOR_ELT(r, i, ScalarReal(s_d));
break;
case CPLXSXP:
for(j=0,in.r=0,in.i=0; j<length(e); j++){
in.r += COMPLEX(e)[j].r;
in.i += COMPLEX(e)[j].i;
}
SET_VECTOR_ELT(r, i, ScalarComplex(in));
break;
default:
SET_VECTOR_ELT(r, i, e);
}
}
setAttrib(r, R_NamesSymbol, getAttrib(lst, R_NamesSymbol));
UNPROTECT(1);
return r;
}これを呼び出す以下のRのコードのファイルexample_list.Rを作ります。
dll <- paste("example_list", .Platform$dynlib.ext, sep = "")
dyn.load(dll)
lst <- list(
i = as.integer(1:5),
r = rnorm(5),
c = c(1 - 2i, 2 + 3i/2),
t = c("a", "b", "c"),
l = list(a=1, b=2, c=3)
)
print("input:")
print(lst)
print("output:")
print(.Call("listsum", lst))
dyn.unload(dll)コンパイルして実行したら、数値データだけ合計されているのが分かります。
R CMD SHLIB example_list.R
R --slave -f example_list.R[1] "input:"
$i
[1] 1 2 3 4 5
$r
[1] -0.8556567 -2.7258396 0.4328710 1.5457363 -0.6162747
$c
[1] 1-2.0i 2+1.5i
$t
[1] "a" "b" "c"
$l
$l$a
[1] 1
$l$b
[1] 2
$l$c
[1] 3
[1] "output:"
$i
[1] 15
$r
[1] -2.219164
$c
[1] 3-0.5i
$t
[1] "a" "b" "c"
$l
$l$a
[1] 1
$l$b
[1] 2
$l$c
[1] 34.3.1
Scalar*関数
Cのソースコード中にも説明を入れましたが、Rにはスカラー型が無く、要素1のベクターの唯一の要素をスカラー変数として扱います。このため、スカラー変数をつくるときすら、まずは型を指定して要素1のベクターをつくり、その唯一の要素に値を代入する必要があって煩雑でした。今は、以下の関数が用意されており、手軽にスカラー変数を作ることが出来ます。
SEXP ScalarReal(double);
SEXP ScalarInteger(int);
SEXP ScalarLogical(int)
SEXP ScalarRaw(Rbyte);
SEXP ScalarComplex(Rcomplex);
SEXP ScalarString(SEXP);
SEXP mkString(const char *);4.3.2 名前でリストの要素を参照する
ベクトルでも可能ですが、リストは名前で要素を参照することが多く、Writing R Extensionsにその典型的なコードが紹介されています。
SEXP getListElement(SEXP list, const char *str)
{
SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);
for (int i = 0; i < length(list); i++)
if(strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {
elmt = VECTOR_ELT(list, i);
break;
}
return elmt;
}これをこのままコピペして、
と言うようなことが出来るわけですが、この関数、標準で提供されないのでしょうか。
4.3.3 Listと言う言葉が指し示すモノ
現在のRでListとされているモノは、昔のRでListとされていたモノとは違い、Cから呼ぶコードにはその痕跡が残っています。上述の例では、isNewListで型チェックをしていた他、allocListではなくallocVectorでListを作成しています。
4.4 表現式を引数にとって、実行する関数
Cのコード中で表現式の評価もできます。
以下の内容のファイルexample_eval.cを作ります。
#include <R.h>
#include <Rinternals.h>
SEXP example_eval(SEXP expr, SEXP rho)
{
SEXP x;
/* 型チェック */
if(!isEnvironment(rho))
error("'rho' should be an environment");
x = PROTECT(allocVector(REALSXP, 3));
REAL(x)[0] = 1.0;
REAL(x)[1] = 1.5;
REAL(x)[2] = 3.1;
/*
環境変数rhoに変数xを定義する
installは文字列をSEXP構造体に変換し、変数xがsymbol tableにないときには付け加える
*/
defineVar(install("x"), x, rho);
UNPROTECT(1);
return eval(expr, rho);
}今回はexample_eval.cのコンパイルをしてから、Rを起動して、
R CMD SHLIB example_eval.c
R --slave -e "system(file.path(R.home(\"bin\"), \"Rgui.exe\"))"以下のように命令して実行してみましょう。
dll <- paste("example_eval", .Platform$dynlib.ext, sep="")
dyn.load(dll)
# エラーが出たときにもDLLがunloadされるようにすると、Rを起動したままデバッグしやすい
tryCatch({
print(.Call("example_eval", quote(sum(x)), new.env()))
}, finally = {
dyn.unload(dll)
})xにc(1.0, 1.5, 3.1)が代入されて、sum(x)が評価されます。
[1] 5.64.4.1
.GlobalEnvを含む任意の環境にオブジェクトをつくれる
環境を渡しているので自明ですが、new.env()を.GlobalEnvに書き換えると、グローバル領域に変数xをつくって計算します。確認しておきましょう。
dll <- paste("example_eval", .Platform$dynlib.ext, sep="")
dyn.load(dll)
# 変数xがあれば消す
if(exists("x")) rm(x)
print(.Call("example_eval", quote(sum(x)), .GlobalEnv))
# 変数xを表示
if(!exists("x")){
print("'x' doesn't exist.")
} else {
print(x)
}
dyn.unload(dll)[1] 1.0 1.5 3.14.5 関数を引数にとって、実行する関数
関数を引数にとることもできます。
これまでと同様にexample_call.cとexample_call.Rを用意します。
example_call.c:
#include <R.h>
#include <Rinternals.h>
SEXP example_call(SEXP x, SEXP fn, SEXP rho)
{
SEXP R_fcall, r;
if(!isNumeric(x))
error("'x' must be numeric.");
if(!isFunction(fn))
error("'fn' must be a function");
if(!isEnvironment(rho))
error("'rho' should be an environment");
/* fnのAbstract Syntax Treeを構築 */
R_fcall = PROTECT(lang2(fn, R_NilValue));
/* 引数にxをセット */
SETCADR(R_fcall, x);
/* ASTを環境rhoで評価する */
r = eval(R_fcall, rho);
UNPROTECT(1);
return r;
}CADRやnilといった語が見えて、Rが関数型プログラミング言語であることが分かります。
example_call.R:
dll <- paste("example_call", .Platform$dynlib.ext, sep="")
dyn.load(dll)
print(
.Call("example_call", as.numeric(1:10), sum, new.env())
)
dyn.unload(dll)example_call.cをコンパイルしてexample_call.Rを実行しましょう。
R CMD SHLIB example_call.c
R --slave -f example_call.R計算できました。
[1] 554.6 関数内部でRの式を組み立てて実行する関数
引数に無くてもRの関数を呼ぶことができます。
これまでと同様にexample_llcall.cとexample_llcall.Rを用意します。
example_llcall.c:
#include <R.h>
#include <Rinternals.h>
SEXP example_llcall(SEXP x, SEXP rho)
{
SEXP R_fcall, ans, s, t;
if(!isNumeric(x))
error("'x' must be numeric.");
if(!isEnvironment(rho))
error("'rho' should be an environment");
/* mean(x, na.rm=TRUE)のAbstract Syntax Treeを構築 */
/* 命令用の領域を確保; tとsには同じポインタが入る */
t = s = PROTECT(allocList(3));
/* LANGSXPを指定 */
SET_TYPEOF(s, LANGSXP);
/* 最初の領域には "mean" を入れる */
SETCAR(t, install("mean"));
/* tを次の領域に移動させる; sは最初の領域のまま */
t = CDR(t);
/* 次の領域には x（へのポインタ）を入れる */
SETCAR(t, x);
t = CDR(t);
/* 最後の領域には論理型でTRUE（非0）を入れて、
名前"na.rm"をつける */
SETCAR(t, ScalarLogical(1));
SET_TAG(t, install("na.rm"));
ans = eval(s, rho);
UNPROTECT(1);
return ans;
}example_llcall.R:
dll <- paste("example_llcall", .Platform$dynlib.ext, sep="")
dyn.load(dll)
x <- 1:11
x[5] <- NA
print(
.Call("example_llcall", x, new.env())
)
dyn.unload(dll)example_llcall.cをコンパイルしてexample_llcall.Rを実行しましょう。
R CMD SHLIB example_llcall.c
R --slave -f example_llcall.R計算できました。
[1] 6.15 DLLで引数の型や呼び出し名を指定する
利用上必須ではないですし大きなメリットも無いのですが、DLLで.Cと.Callと.Fortranと.Externalで呼び出し名をつけることができ、.Cと.Fortranに型制約をつけられます7。
試しに関数repetitivenessをtautologyに、minusplusをplusminusに置き換えてみましょう。
example_nr.c:
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Rdynload.h>
static double minusplus(double *x, int *len)
{
int i, z=1;
double y = 0;
for(i=0; i<*len; i++){
y = z * x[i];
z = -1 * z;
}
return y;
}
static SEXP repetitiveness(SEXP obj){
return obj;
}
static R_NativePrimitiveArgType minusplus_type[] = {
REALSXP, INTSXP
};
static const R_CMethodDef cMethods[] = {
{"plusminus", (DL_FUNC) &minusplus, 2, minusplus_type},
NULL
};
static const R_CallMethodDef callMethods[] = {
{"tautology", (DL_FUNC) &repetitiveness, 1},
NULL
};
/* `example_nr`はsuffixなしファイル名の部分なので、ファイル名が変われば変わる */
void R_init_example_nr(DllInfo *info)
{
/* .Fortran，.Externalの定義があれば、第4，第5引数に指定する */
R_registerRoutines(info, cMethods, callMethods, NULL, NULL);
}example_nr.R:
dll <- paste("example_nr", .Platform$dynlib.ext, sep="")
a <- dyn.load(dll)
print(
.Call("tautology", 443)
)
x <- 1:10
print(
.C("plusminus", as.numeric(x), as.integer(length(x)))
)
dyn.unload(dll)これまで通り、example_nr.cをコンパイルして、example_nr.Rを実行します。
R CMD SHLIB example_nr.c
R --slave -f example_nr.R計算できました。
[1] 443
[[1]]
[1] 1 2 3 4 5 6 7 8 9 10
[[2]]
[1] 10なお、不正な型を入れると、
型があわないとエラーが出て止まります。
.C("plusminus", as.numeric(x), 3.1) でエラー:
2 引数に対する不正な型です（plusminus への呼び出しにおいて）
呼び出し: print -> .C
実行が停止されました6
可変長引数のときに使う.External
.Callでだいたい間に合うと思いますが、.Callは可変長引数...をとることができません8。.Externalは可変長引数を取れる呼び出しで、Cの関数ではCARで可変長引数の最初の要素を参照し、CDRで可変長引数の最初の要素を消すことを、可変長引数がR_NilValueになるまで繰り返すLispの中身感のあるコーディングを行ないます。具体的例はWriting R
ExtensionsのSEXP showArgs(SEXP args)関数が過不足ないので、参照してください。
7 C言語から利用できるRのAPI
RからCの関数を呼び、RとCでデータをやり取りし、CでRのデータ型を操作するための機能を提供するだけではなく、Rで標準で使える積分と最適化を含む数学・確率統計関数やソートなどの関数が提供されています。Writing R Extensionsの[6 The R API: entry points for C code]に列挙されているので目を通されるとよいと思いますが、利用頻度の高そうなものを紹介します。
7.1
Rコンソール出力関数REprintf
REprintfはC言語のprintfと同様に使える、Rコンソール出力関数です。Windows版Rで開発するとgdbの利用が制限されるので、printfデバッグforeverと重宝するかもです。
7.2 メモリー管理関数群
R_allocは一時メモリーを割り当ててもらう関数で、エラーの発生やユーザーによる中断を含めて、.Cと.Callと.ExternalによるC呼び出しが終了したら、自動的にメモリーが開放され8ます。メモリーリークの防止に有用です。
R_Calloc/R_Realloc/R_FreeはRを通してメモリーを割り当ててもらう関数です。昔のドキュメントではCalloc/Freeが説明されていたのですが非推奨になっており、#define STRICT_R_HEADERSをしてから使うとエラーになります。GitHubにあげているあれが･･･。
これらの関数は、double *a = (double *)R_alloc(100, sizeof(double));と言う風に、対応するC言語の標準関数と同様に使えます9。
7.3 ユーザー中断の許可
Cのコードの実行中はユーザー中断が許可されないのですが、ループ中にR_CheckUserInterrupt(void);を挟んでおくと、CTRL+Cなどが押されていれば、そこで中断されます。
7.4 乱数生成関数
マルコフ連鎖モンテカルロ法などで必須になる乱数です。C言語のメルセンヌ・ツイスタのライブラリを探してくるのは難しくないわけですが、Rに内蔵されているものを使うこともできます。
/* 乱数利用開始前に必ず */
void GetRNGstate(void);
/* user_unif_init(Int32)に変わるseed設定関数 */
double R_unif_index(double);
/* 正規化された一様分布，正規分布，指数分布からの乱数 */
double unif_rand();
double norm_rand();
double exp_rand();
/* 乱数利用終了時に必ず唱える（seed更新） */
void PutRNGstate(void);R_ext/Random.hに定義されているので、使う前にはincludeしてください。なお、dnorm/pnorm/qnorm/rnormと言った確率統計に関するRの同名の関数とほぼ同様のC関数も提供されていて、Rmath.hに定義されています。
7.5 数学関数/定数
Rmath.hを見ると、実装はもちろんライブラリを探すのも手がかかる最適化と積分も含めて、Rの基本パッケージで使える関数と定数と同様の関数と定数がCに提供されていることが分かります。
7.5.1 最適化
optim関数の中にあるNelder Mead法、BFGS法、Conjugate
gradients法、L-BFGS法、Simulated
annealing法それぞれの関数nmmi, vmmin,
cgmin, lbfgsb,
saminを呼ぶことができ、最適化問題を解けます。これらの手法につきもののエラーや警告の扱いはR_tryCatchError,
R_tryCatch,
R_withCallingErrorHandler関数で行なえます。
試しにL-BFGS法で、正規分布を仮定して平均と標準偏差を推定してみましょう。
example_lbfgs.c:
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Applic.h>
#include <math.h>
/* データセットを表す構造体 */
typedef struct {
unsigned int n;
double *y;
} DATASET;
double objf(int n, double *par, void *ex)
{
/* 正規分布の対数尤度関数（なのでnは2で決まっているので使わない） */
DATASET *ds = ex;
double theta = par[0];
double mu = par[1];
double s, d;
int i;
for(s=0,i=0; i<ds->n; i++){
d = ds->y[i] - mu;
s += d*d;
}
return -1.0*(-1.0*(double)(ds->n)*log(theta*theta)/2 - s/(theta*theta)/2);
}
void objfg(int n, double *par, double *gr, void *ex)
{
/* objfのグラディエントの計算; 一般的には数値微分した方が使い勝手は良い */
DATASET *ds = ex;
double theta = par[0];
double mu = par[1];
double s, d, res;
int i;
for(s=0,res=0,i=0; i<ds->n; i++){
d = ds->y[i] - mu;
res += 2*d;
s += d*d;
}
gr[0] = (double)(ds->n)/theta - s/(theta*theta*theta);
gr[1] = -1.0*res/(theta*theta)/2;
}
SEXP mle_nd(SEXP init, SEXP y)
{
DATASET ds;
int n, i;
SEXP r;
double *x, lower=-1e+10, upper=1e+10, Fmin;
int fail=0, fncount=0, grcount=0, *nbd;
char *msg = R_alloc(60, sizeof(char));
if(!isVector(init))
error("'x' should be a vector.");
if(!isReal(init))
error("'x' should be of numeric.");
if(!isVector(y))
error("'y' should be a vector.");
if(!isReal(y))
error("'y' should be of numeric.");
n = length(init);
x = (double *)R_alloc(n, sizeof(double));
nbd = (int *)R_alloc(n, sizeof(int));
for(i=0;i<n;i++){
x[i] = REAL(init)[i]; /* set the initial value */
nbd[i] = 0; /* unbounded */
}
n = length(y);
ds.n = n;
ds.y = (double*)R_alloc(n, sizeof(double));
for(i=0;i<n;i++){
ds.y[i] = REAL(y)[i];
}
lbfgsb(length(init), 5, x, &lower, &upper, nbd, &Fmin, objf, objfg,
&fail, &ds, 1e+7, 0, &fncount, &grcount, 100, msg, 0, 10);
PROTECT(r = allocVector(REALSXP, 2));
REAL(r)[0] = x[0]; /* sd */
REAL(r)[1] = x[1]; /* mean */
UNPROTECT(1);
return r;
}Writing R
Extensionsのアドバイス通りに、パラメーターはoptimのヘルプを見ながらRの省略値にあわせつつ、nbdとmsgはソースコードにあったコメントの記述を参考に入れました。
example_lbfgs.R:
dll <- paste("example_lbfgs", .Platform$dynlib.ext, sep = "")
dyn.load(dll)
tryCatch({
y <- rnorm(30, mean=1, sd=2)
init <- c(2, 2)
print(.Call("mle_nd", init, y));
}, finally = {
dyn.unload(dll)
})これまで通り、example_lbfgs.cをコンパイルして、example_lbfgs.Rを実行します。
R CMD SHLIB example_lbfgs.c
R --slave -f example_lbfgs.R計算できました。
[1] 2.1553354 0.9978105なお、乱数でデータセットを作っているので、推定量は動かす度に変わります。
7.5.2 積分関数
Gauss–Legendre法のような積分アルゴリズムを考えたくないときには10、Rのintegrate関数の中にある関数で積分することができます。
example_integrate.c:
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Applic.h>
/*
n個の点における被積分関数の値を引数のポインターに書き込むユーザー定義函数
*/
void vectorized_integrand(double *x, int n, void *ex)
{
int i;
for(i=0; i<n; i++){
/* 区間[0 1)で11.21をとる被積分関数を例にしている */
if(x[i]<0)
x[i] = 0;
else if(x[i]>=1)
x[i] = 0;
else
x[i] = 11.21;
}
}
SEXP example_integral(void)
{
/* 入力 */
void *ex = NULL; /* 定義する被積分関数が内部で使うユーザー定義の補助データの構造体へのポインター */
double a = -2, b = 2; /* 区間[-2 2)を積分 */
double epsabs, epsrel; /* 絶対/相対精度 */
/* 作業領域とその大きさ */
int limit, lenw, *iwork;
double *work;
/* 出力 */
int last; /* the output components value, abs.err and subdivisions of the R function integrate */
int neval, ier; /* 不明, 被積分関数が評価された回数, エラーコード */
double result = NAN, abserr = NAN; /* 積分値, 絶対誤差 */
/* Rのintegrate関数にあわせて入力値を設定 */
epsabs = epsrel = sqrt(sqrt(DBL_EPSILON));
limit = 100; /* limitはRのintegrate関数のsubdivisionsに該当*/
lenw = 4 * limit;
iwork = (int *) R_alloc(limit, sizeof(int));
work = (double *) R_alloc(lenw, sizeof(double));
/* 有界区間の積分関数 */
Rdqags(vectorized_integrand, ex,
&a, &b,
&epsabs, &epsrel,
&result, &abserr, &neval,
&ier,
&limit, &lenw, &last,
iwork, work);
return ScalarReal(result);
}連続で積分するときのオーバーヘッドを抑えるためか、作業領域が使いまわせる仕様です。なお非有界区間の場合はRdqagi関数を用います。
example_integrate.R:
dll <- paste("example_integral", .Platform$dynlib.ext, sep = "")
dyn.load(dll)
tryCatch({
print(.Call("example_integral"))
}, finally = {
dyn.unload(dll)
})これまでと同じようにexample_integrate.cとexample_integrate.Rを作成したディレクトリで、
R CMD SHLIB example_integrate.c
R --slave -f example_integrate.Rとすると、自明な積分値が出力されます。
[1] 11.217.6 その他
ソートや文字コード変換、一時ファイル名の取得のAPIも提供されています。
8 OpenMPによる並列化
何重にもループするような処理を高速化、さらに並列処理をする例をあげます。
.Callではなくて.Cを使っているのは、Rのオブジェクトがスレッドセーフで無いためコードがさらに冗長になるためです。
OpenMPと聞くと高度な技能が要求される気がするかも知れませんが、難しく煩雑なことはコンパイラがやってくれるので利用者は楽ができる技術です。
以下の例ではforの前に1行、OpenMPでループを並列化すること、変数a,
b,
mは並列スレッド間で同期処理を行なわないことを宣言しているだけです。
example_openmp.c:
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Rdynload.h>
void gcd(int *len, int *v, int *w, int *r)
{
int i, a, b, m;
#pragma omp parallel for private(a, b, m)
for(i=0; i<*len; i++){
a = v[i];
b = w[i];
while(0!= (m = a % b)){
a = b;
b = m;
}
r[i] = b;
}
}
/* 以下は無くても型チェックが弱くなるだけ */
static R_NativePrimitiveArgType gcd_type[] = {
INTSXP, INTSXP, INTSXP, INTSXP
};
static const R_CMethodDef cMethods[] = {
{"gcd", (DL_FUNC) &gcd, 4, gcd_type},
NULL
};
void R_init_example_openmp(DllInfo *info)
{
R_registerRoutines(info, cMethods, NULL, NULL, NULL);
}C/OpenMPで書いた最小公倍数の計算コードと、Rで書いたものの処理時間を比較します。
example_openmp.R:
dll <- paste("example_openmp", .Platform$dynlib.ext, sep="")
dyn.load(dll)
gcd_c <- function(v, w){
if(length(v)!=length(w)) stop("v and w need to have the same length.")
r <- .C("gcd", length(v), as.integer(v), as.integer(w), integer(length(v)))
r[[4]]
}
gcd_r <- function(v, w){
if(length(v)!=length(w)) stop("v and w need to have the same length.")
r <- integer(length(v))
for(i in 1:length(v)){
a <- v[i]
b <- w[i]
while(0 != (m = a %% b)){
a <- b
b <- m
}
r[i] <- b
}
r
}
n <- 10^7
v <- round(runif(n, min=0.5, max=n-0.5))
w <- round(runif(n, min=0.5, max=n-0.5))
tryCatch({
print(".C and OpenMP:")
print(system.time({
gcd_c(v, w)
}))
print("R:")
print(system.time({
gcd_r(v, w)
}))
}, finally = {
dyn.unload(dll)
})これまでと同じようにexample_openmp.cとexample_openmp.Rを作成したディレクトリで、
R CMD SHLIB example_openmp.c
R --slave -f example_openmp.RとするとDLLの作成と実行ができます。準備が不十分でMakevars/Makevars.winが適切に作られていない場合、OpenMPが有効にならずにコンパイルエラーが出るので注意してください。
[1] ".C and OpenMP:"
ユーザ システム 経過
0.18 0.02 0.07
[1] "R:"
ユーザ システム 経過
5.72 0.00 5.84Rより.Cのコードの方が圧倒的に速いですね。また、.Cはユーザ時間よりも経過時間の方が小さいので、OpenMPによる平行処理が出来ていることがわかります。
ユーザ時間が経過時間よりも短い場合、並列処理はされていません。環境変数OMP_NUM_THREADSに2以上の値を設定してから、再度、実行してみてください。
9 Rのゾンビがメモリに残ってファイルを放さない！
Rから呼んでいるCのコードが異常終了するとRも心中することになるのですが、そのときRのゾンビプロセスが残って、メモリーが圧迫される他、RやCのコードが開いていたファイルやフォルダーがロックされて更新や削除ができなくなったりします。
複数のゾンビがあって始末するのが手間なときは、Linuxであればpkill -9 Rをシェルで、Windowsであればtaskkill /IM Rterm.exe /Fをコマンドプロンプトで入力すればまとめて始末できます。
10 まとめ
Rからデータを貰って、Rに処理結果を返すだけなので、Rから呼ぶCのコードは（Rに新たなUIを提供する目的でなければ）UIを書かなくて済みますし、さらにライブラリが色々と提供されているので、CでOSのAPIを叩いてアプリケーションを書くよりは楽です。
そして有用さが確認できたら（本稿では説明しませんでしたが）すかさずパッケージ化できます。もはやRはC言語の開発フレームワークの一つと見做せます。
数値解析の高速化のためと考えるとRcppの方が労力対効果で効率的ですが、ライブラリ的なものを整備するのには.Callと.Externlの方がRcppが不要になるので便利な場面もあります。
Rを使っているCプログラマであれば、試して見る価値は大きいです。
1990年代中盤から広範に一世を風靡したPerlだってXS言語と言うC言語を埋め込む仕込むがありましたし、Fortranも2003からCから呼ばれたり、Cを呼ぶための機能拡張がされました。Python，Ruby, Java, JuliaからもC-callは可能です。↩︎
定型的な決まりきった記述。↩︎
線形代数の演算子が使えること、最適化をしてくれたり疎行列の固有値を戻す線形代数ライブラリが使える点が大きな違いになります。↩︎
2001年1月版の翻訳もあるのですが、例示されているコードがまだ動く一方で、21年の間に加筆された部分もそれなりあり、古くなった記述もあります。↩︎
%R_HOME%は標準インストール時はC:\Program Files\R\R-4.2.2\binになります。↩︎coerceVectorでSEXP型の内部型の変換時に、PROTECT_INDEX型の変数を用意しておき、PROTECT_WITH_INDEXマクロでインデックスを保存、REPROTECTマクロでインデックスを指定してスタックの中身を入れ替える例が示されていました。↩︎Fortranに関してはヘッダーファイルを見ただけで、実際のコードで確認していません。
externしたFortranのサブルーチンを指定することになるのだと思いますが。↩︎Rの関数で可変長引数をリストに変換できるので、困ることは少ないと思いますが。↩︎
厳密にはちょっと違っていて、なぜか
void *ではなくchar *で定義されているので型キャストがいります。↩︎パフォーマンス向上にはより良い積分アルゴリズムを探してくるのが重要ですが、限界を攻めないときはこれで十分だと思います。↩︎