이 글에 사용된 코드와 명령은 MacOSX와 Linux에서 확인된 것이다. 특별히 설치되어야 하는 것은 없다. 기본적으로 R과 기본 개발도구가 설치되어 있으면 된다. MacOSX에서는 Xcode를 설치하면 된다. Fortran/C/C++ 무엇으로 짜야 하나 고민이 된다. 다음은 확인된 사실은 아니지만...
- R에서 vectorize할수 있는 데까지 다 했는데도 속도를 빠르게 하고 싶다면 Fortran/C/C++을 생각해야 한다.
- Fortran이 C보다 우월한가? 논쟁만 부르고 결론은 나지 않는 질문이다. 그럴 수도 있고 아닐 수도 있다.
- R에서는 공식적(?)으로 Fortran 77만 지원한다. Fortran 95로 작성한 코드를 R에서 부르는 방법이 없는 건 아니지만... 어떤 사람은 Fortran 77은 절대로 쓰지 말아야한다고도 한다. Fortran 95도 있고 최근 2010년에 표준화된 흔히 부르는 이름으로 "Fortran 2008"도 있다.
.Call
R 패키지를 개발할 때나 아니면 그냥 함수를 짤 때 C/C++/Fortran으로 짠 코드를 불러서 쓸 수 있다. R 함수 중.C()는 컴파일된 C 코드를 호출할 때, .Fortran()는 컴파일된 Fortran 77 코드를 호출할 때 쓸 수 있다. 이것보다 좀더 유연하고 성능도 나은 것이 .Call()이다. 이것은 C/C++ 코드에 R 객체를 넘겨줄 때 사용된다. .External()도 마찬가지인데 .Call()은 고정된 수의 인자를 넘기고 .External()은 argument list로 넘긴다는 데에 차이가 있다.
.Call()을 통해 사용할 C 코드를 작성할 때-
R.h,Rinternals.h헤더를 반드시 include한다. - R 객체를 받을 때
SEXP,REALSXP,INTSXP,STRSXP타입을 이용한다. - R 객체의 속성은
getAttrib()로 뽑아낼 수 있다. - R 객체 안에 들어있는 데이터에 대한 포인터는
REAL(),INTEGER()등을 이용하여 얻을 수 있다. - R 객체를 만들 때는
allocVector(),allocMatrix()를 이용하여SEXP타입 객체로 만든다. - 이때 항상
PROTECT()와UNPROTECT()를 이용한다. 안전한 메모리 확보를 위한 장치이다.
예를 들어 다음과 같은 함수를 작성하고
myfun.c 파일에 저장하였다고 하자.#include <R.h>
#include <Rinternals.h>
SEXP say(SEXP x) {
return x
}
컴파일은 다음과 같이 한다.
$ R CMD SHLIB myfun.c컴파일하면 만들어지는
myfun.so를 R에서 불러 사용한다. R에서 다음과 같이 로딩한다.> dyn.load("myfun.so")
로딩이 되었으면 다음처럼 실행할 수 있다.> .Call("say", 3)
[1] 3
예제: 평균 계산하기
이제 C로 평균을 구하는 함수를 다음과 같이 작성하고mystat.c 파일로 저장하였다고 하자.// mystat.c
#include <R.h>
#include <Rinternals.h>
SEXP mean(SEXP x) {
double *vec = REAL(x);
int n = length(x);
SEXP val = PROTECT(allocVector(REALSXP, 1));
double sum = 0;
for(int i=0; i < n; i++) sum += vec[i];
REAL(val)[0] = sum / n;
UNPROTECT(1);
return val;
}
컴파일을 하면 mystat.so 파일이 생긴다.
$ R CMD SHLIB mystat.cR에서 이것을 불러쓰기 위해 다음처럼 wrapper를 작성하자.
mymean <- function(x) {
if(!is.loaded("mystat")) dyn.load("mystat.so")
.Call("mean", as.numeric(x))
}
이제 R에서 mymean()이라는 함수를 사용하면 된다.
속도 비교
R에는 물론 평균을 계산해주는mean() 함수가 있다. 만약 R에서 for으로 평균을 구하면 어떨까? iMac 3.06GHz Intel Core 2 Duo에서 1e7개의 난수를 더한 후 평균을 구하는 데에 다음처럼 오랜 시간이 걸린다.
> n <- 1e7
> x <- runif(n)
> system.time( {
+ total <- 0
+ for(i in 1:n) total <- total + x[i]
+ total / n
+ })
user system elapsed
8.436 0.076 9.165
물론 R의 mean() 함수를 이용하면 빠르게 계산된다.
> system.time( mean(x) ) user system elapsed 0.079 0.001 0.086방금 C로 짠 함수의 경우에는:
> system.time( mymean(x) ) user system elapsed 0.020 0.000 0.023R의 내장
mean()보다 빠른 속도를 보여주고 있는데 데이터에 따라서 다를 수 있으니 항상 성능이 좋다고는 못한다. 어쨌든 R에서 for문 돌리는 것에 비하면 빠르다.
행렬
R에서 넘어온 행렬은 컬럼을 쌓은 벡터로 받아서 사용한다. $p\times q$ 행렬 $A$와 $q\times r$ 행렬 $B$의 행렬곱은 다음처럼 계산할 수 있다.#include <R.h>
#include <Rinternals.h>
/**
* @SmatA Sp by Sq matrix
* @SmatB Sq by Sr matrix
**/
SEXP prod(SEXP SmatA, SEXP SmatB, SEXP Sp, SEXP Sq, SEXP Sr) {
int p = INTEGER(Sp)[0];
int q = INTEGER(Sq)[0];
int r = INTEGER(Sr)[0];
// R matrix as a vector stacked column by column
double* a = REAL(SmatA);
double* b = REAL(SmatB);
SEXP ans = PROTECT(allocVector(REALSXP, p*r));
double sum;
int i, j, k;
for(i=0; i < p; i++) {
for(k=0; k < r; k++) {
sum = 0;
for(j=0; j < q; j++) {
sum += a[i + j*p] * b[j + k*q];
}
REAL(ans)[i + k*p] = sum;
}
}
UNPROTECT(1);
return ans;
}
