俺俺lapply&sapplyを書いてみる【Rの中身はlisp編】

というかWriting R Extensionsにある例をやってみているという感じなんですけど。

C言語が(も)あんまり分かっていない僕ですが、それでもRのベクトルや行列のオブジェクトをCレベルで触るにはどうすればいいか段々分かってきました。こうなってくるとRの色々なオブジェクトをCレベルで見ていきたくなります。で、触りたくなったのが「関数オブジェクト」。applyとかsapplyの中でどういう風に持たれているかが知りたくなったんですね。

というわけでコードを書いてみた。もうなんというかcadrとか出てきて完全にlispな乗りです。

#include <R.h> 
#include <Rinternals.h> 

SEXP mylapply(SEXP list, SEXP fn, SEXP rho) 
{ 
  // evalとかの関係で環境も引数に渡してあげないといけない
  int i, n = length(list); 
  SEXP R_fcall, ans; 
  if(!isNewList(list)) error("‘list’ must be a list"); 
  if(!isFunction(fn)) error("‘fn’ must be a function"); 
  if(!isEnvironment(rho)) error("‘rho’ should be an environment"); 
  
  PROTECT(R_fcall = lang2(fn, R_NilValue)); // 2つのSEXPをconsして、そのTYPEOFをLANGSXPにする
  Rprintf("Value of fn : ");
  PrintValue(fn);
  PROTECT(ans = allocVector(VECSXP, n)); 
  for(i = 0; i < n; i++) { 
    Rprintf("i = %d\n",i);
    Rprintf("Before : R_fcall -> "); // 評価はまだされていない
    PrintValue(R_fcall);
    SETCADR(R_fcall, VECTOR_ELT(list, i)); // S式の肝。(+ 1 2)というような形のリストにする
    Rprintf("After : R_fcall -> ");
    PrintValue(R_fcall);
    Rprintf("CAR -> ");
    PrintValue(CAR(R_fcall));    
    Rprintf("CDR -> ");
    PrintValue(CDR(R_fcall));
    Rprintf("CAR(CDR) -> ");
    PrintValue(CAR(CDR(R_fcall)));
    SET_VECTOR_ELT(ans, i, eval(R_fcall, rho)); // ans[i]にeval(R_fcall, rho)の結果を格納
  } 
  Rprintf("Attribution : ");
  PrintValue(getAttrib(list, R_NamesSymbol)); // 名前がついているlistにする
  setAttrib(ans, R_NamesSymbol, getAttrib(list, R_NamesSymbol)); 
  UNPROTECT(2); 
  return(ans); 
} 

SEXP mysapply(SEXP vec, SEXP fn,  SEXP rho)
{
  // mylapplyを使えばいいんだけど、わざと書いてみる
  int i, n = length(vec);
  SEXP R_fcall, ans;
  PROTECT(R_fcall = lang2(fn, R_NilValue));
  PROTECT(ans = allocVector(VECSXP, n));
  // VECTOR_ELT(vec, i);とはできない(vectorだから)。listじゃないとだめなようだ。
  SETCADR(R_fcall, vec);
  PrintValue(R_fcall);
  ans = eval(R_fcall, rho);
  UNPROTECT(2); 
  return(ans); 
}

このCのコードをいつものごとく

R CMD SHLIB fuga.c

とコンパイルしてやる。で、以下のコードで実行してやるんだけど、昨日のところのシンボルオブジェクトみたいなのが出てきてどう考えてもlispです。

(+ 1 2 3)

こんな感じ。

> a <- list(a = 1:5, b = rnorm(3))
> dyn.load("/tmp/fuga.so")
> .Call("mylapply", a, sum, new.env()) 
Value of fn : function (..., na.rm = FALSE)  .Primitive("sum")
i = 0
Before : R_fcall -> .Primitive("sum")(NULL)
After : R_fcall -> .Primitive("sum")(1:5)
CAR -> function (..., na.rm = FALSE)  .Primitive("sum")
CDR -> [[1]]
[1] 1 2 3 4 5

CAR(CDR) -> [1] 1 2 3 4 5
i = 1
Before : R_fcall -> .Primitive("sum")(1:5)
After : R_fcall -> .Primitive("sum")(c(-0.686190593452094, -0.709346104248138, -0.56883384303222
))
CAR -> function (..., na.rm = FALSE)  .Primitive("sum")
CDR -> [[1]]
[1] -0.6861906 -0.7093461 -0.5688338

CAR(CDR) -> [1] -0.6861906 -0.7093461 -0.5688338
Attribution : [1] "a" "b"
$a
[1] 15

$b
[1] -1.964371

> .Call("mysapply", 1:5,function(x){x^3}, new.env()) 
function(x){x^3}
(1:5)
[1]   1   8  27  64 125

Rinlinedfuns.h

lispな正体を見ていく。下のコードはconsとかやっていってTYPEOFをLANGSXPとかにしている。lang4とかもどんどんconsしているだけ。全然大したことない。

INLINE_FUN SEXP lcons(SEXP car, SEXP cdr)
{
    SEXP e = cons(car, cdr);
    SET_TYPEOF(e, LANGSXP);
    return e;
}

INLINE_FUN SEXP lang1(SEXP s)
{
    return LCONS(s, R_NilValue);
}

INLINE_FUN SEXP lang2(SEXP s, SEXP t)
{
    PROTECT(s);
    s = LCONS(s, list1(t));
    UNPROTECT(1);
    return s;
}

INLINE_FUN SEXP lang3(SEXP s, SEXP t, SEXP u)
{
    PROTECT(s);
    s = LCONS(s, list2(t, u));
    UNPROTECT(1);
    return s;
}

INLINE_FUN SEXP lang4(SEXP s, SEXP t, SEXP u, SEXP v)
{
    PROTECT(s);
    s = LCONS(s, list3(t, u, v));
    UNPROTECT(1);
    return s;
}

listをただただconsしていくlist1、list2、list3、list4というのも用意してある。

/* Shorthands for creating small lists */

INLINE_FUN SEXP list1(SEXP s)
{
    return CONS(s, R_NilValue);
}


INLINE_FUN SEXP list2(SEXP s, SEXP t)
{
    PROTECT(s);
    s = CONS(s, list1(t));
    UNPROTECT(1);
    return s;
}


INLINE_FUN SEXP list3(SEXP s, SEXP t, SEXP u)
{
    PROTECT(s);
    s = CONS(s, list2(t, u));
    UNPROTECT(1);
    return s;
}


INLINE_FUN SEXP list4(SEXP s, SEXP t, SEXP u, SEXP v)
{
    PROTECT(s);
    s = CONS(s, list3(t, u, v));
    UNPROTECT(1);
    return s;
}

新しいガベージコレクション機能

Original Replacement
foo = VECTOR(bar)[i] foo = VECTOR_ELT(bar, i)
VECTOR(foo)[j] = bar SET_VECTOR_ELT(foo, j, bar)
foo = STRING(bar)[i] foo = STRING_ELT(bar, i)
STRING(foo)[j] = bar SET_STRING_ELT(foo, j, bar)