R 프로그래밍: R에서 C 함수 부르기

R 프로그래밍: R에서 C 함수 부르기

이 글에 사용된 코드와 명령은 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.c

R에서 이것을 불러쓰기 위해 다음처럼 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.023 

R의 내장 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;
}