본문으로 바로가기



Make your own R package! - 7. R과 C 언어 연결하기 / R CMD SHLIB 명령어 사용 및 dyn.load 함수

category R Programming/Rpackage 2017.08.03 18:25
Make your own R package! - 7. R과 C 언어 연결하기

오늘의 주제는 R와 C 언어를 연결하기 위한 준비단계로 이전 포스팅에서 짰던 dpareto함수를 C언어로 짜는 것을 알아보도록 하겠다. R을 주로 포스팅하는 필자가 C언어에 대하여 포스팅하는 이유는 C언어를 다룰수 있으면 R이 제일 취약한 부분인 루프(loop) 부분을 C로 넘겨서 계산을 할 수 있도록 연결해줄 수 있기 때문이다.

C 언어로 함수 정의하기

이전 포스팅에서 다루었던 파레토 분포의 확률밀도함수 dpareto를 C언어를 사용하여 만들어 보도록 하자. 먼저 가장 기본이 되는 변수 하나씩 받을 수 있는 함수를 paretodens라는 이름으로 다음과 같이 정의할 수 있다. C의 경우 함수의 앞에 함수의 결과값이 어떠한 형식을 띄는지 명시해주기 때문에 paretodens앞에 double이 붙어있으며, 각 입력값과 함수 안에서 정의되는 변수의 형태 역시도 앞에 명시해 주도록 되어있다.

double paretodens( double x, double alpha, double x_m, int log_tag ) {
  double density;
  
  if (alpha <= 0.0 || x_m <= 0```{r cache=TRUE}```{r cache=TRUE}.0){
    printf ("The values of all parameters should be positive.\n");
    density = -1.0;
  } else {
    if ( x >= x_m ) {
      density = log(alpha) + alpha * log(x_m) - (alpha + 1.0) * log(x);
    } else {
      density = log(0.0);
    }
  }
  
  if (log_tag != 1) {
    return density;
  } else {
    return exp(density);
  }
}

위의 paretodens 함수의 경우 처음 모수가 적절한 값을 갖는지에 대하여 체크를 한 후, 적합하지 않는 경우 density값을 음수값으로 설정{r cache=TRUE}{r cache=TRUE}한다. 그 외 모수의 적합도 테스트를 통과한 다음에는 정의역 안에서는 밀도함수 값을 로그 스케일에서 계산을 하고, 그외의 영역에서는 -Inf값을 부여해서 나중에 exp함수를 거쳐 0이 되도록 하는 방식을 짜여져 있다. 이렇게 짠 paretodens 함수를 main 함수를 통해서 시스템의 커널을 통해서 입력값을 전달받아 paretodens 함수로 넘겨주는 함수를 다음과 같이 설정한다.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* declare function */
double paretodens( double x, double alpha, double x_m, int log );
double x, alpha, x_m, value;
char line [100];
int log_tag;

int main(void) {
  /* printf( " Enter values for x , alpha, x_m : " ); */
  fgets( line , sizeof( line ), stdin );
  sscanf( line , "%lf %lf %lf ", &x , &alpha , &x_m );
  value = paretodens(x, alpha, x_m, 1);
  
  if ( value < 0 ) {
    printf( "Error!\n" );
  } else {
    printf( "%.5f\n", value );
  }
  return 0;
}

위의 코드에서 stdio.hstdlib.h, 그리고 math.h는 C에서 사용하는 헤더파일이다. 필자의 경우는 그냥 우리가 R에서 기본 베이스 패키지를 불러와 이미 정의된 함수를 사용하는 것처럼, C에서 역시 미리 정의된 함수들을 불러와서 사용한다고 이해하고 있다. ’stdio.h의 경우는 입출력을 담당하는 함수를 불러오기위하여 추가시킨다.’는 정도로만 이해하고 넘어간다.

우리가 사용할 파일의 이름을 dpareto.c라고 가정하면, 파일은 다음 첨부파일과 같이 main함수로 들어가기전 사용되는 함수의 키워드가 나오고, main함수 그리고 그 뒤로 paretodens가 정의되어 있는 구조가 되어야 한다. 이렇게 생성된 파일을 리눅스 시스템을 기준을 다음의 명령어를 사용하여 컴파일을 한 후, 실행파일로 만들 수 있다.

dpareto.c

gcc -o dpareto dpareto.c -lm
make dpareto

컴파일이 끝난 후 이것을 실행파일로 만들어 실제 함수를 작동시키기 위해서는 make 명령어를 이용해서 실행가능 파일로 만들어줘야한다. 이 과정을 마치면 시스템 명령어를 사용하여 다음과 같이 함수의 작동을 확인 할 수 있다.(아래의 경우는 리눅스 상에서의 명령어이다.)

echo -4.0 2.0 3.0 | ./dpareto
echo  2.0 4.0 3.0 | ./dpareto

.C 를 사용한 R과의 연결

현재 R과 C를 연결하는 방법은 크게 두 가지, .C.Call 방법이 있다. 오늘은 비교적 단순한 방법인 .C를 사용하여 R에서 C언어 이용하는 방법에 대하여 알아보도록하자. .C는 R에서 사용되는 C언어 함수에 다음의 세가진 제약조건을 부여한다.

  1. R에서 사용되는 C함수는 void를 반환한다.
  2. C언어 함수로 전달되는 함수 인자들은 모두 pass by referece 개념을 사용하여 전달한다. 따라서 정의된 함수의 입력값들이 모두 포인터를 사용하여 전달된다.
  3. C언어로 정의된 함수는 모두 <R.h> 헤더 파일을 포함하고 있어야 한다.

위의 세가지 조건을 만족시키기 위해서는 우리가 사용하는 함수는 (예를 들어 paretodens라고 하자.) 다음과 같은 형식을 사용해야만 한다.

#include <R.h>
...
void paretodens(double *x, double *alpha,
                double *x_m, int *mlength,
                double *density){
...
}

함수의 반환값이 void이므로, 결과값을 받아올 변수까지 넘겨주는 형식을 취하게 된다(위의 코드에서는 density가 값을 받아오는 역할을 한다.) 먼저 R코드에서 어떻게 .C를 사용하여 C언어로 짜여진 함수를 불러오는지에 대하여 알아보자. 함수의 입력값에서 mlenght는 결과값을 받을 density 배열의 크기 정보를 담고 있는 변수이다. 함수의 입력값을 포인터로 받기 때문에 paretodens 함수 안에서는 sizeof 함수를 사용하여 변수의 크기를 알아낼 수 없으므로, 다른 값들과 함께 전달해줘야 한다.

R에서 .C를 통해 C언어 함수 불러오기

먼저 어떻게 C로 짜여진 함수를 R에서 불러올 수 있는지에 대하여 알아보자. 우리가 정의한 C언어 함수의 이름이 paretodens라고 할 때, 다음과 같이 .C 함수를 사용하여 paretodens 함수에 값을 넘겨줄 수 있다. 아래의 코드는 이전 포스팅인 함수의 벡터화에서 정의한 코드를 응용하였다.

dpareto <- function(x, alpha, x_m, log = FALSE) {
  mlength <- max(length(x), length(alpha), length(x_m))

  # Input vectorizing
  density <- vector(mode="numeric", length = mlength)
  x     <- rep(x, length.out = mlength)
  alpha <- rep(alpha, length.out = mlength)
  x_m   <- rep(x_m, length.out = mlength)

  # Call .C
  result <- .C("paretodens",
               as.double(x),
               as.double(alpha),
               as.double(x_m),
               as.integer(mlength),
               value = as.double(density))$value
  
  if (any(is.nan(result))) {
    warning("All parameters should be positive.\n NaN is generated.\n")
  }
  
  if (log) result else exp(result)
}

.C를 통하여 전달되는 입력값들이 모두 길이가 같은 벡터라는 것에 주목하자. 그리고 paretodens를 통하여 계산을 한 후의 결과값을 R은 리스트 형태로 받게 되는 것을 알 수 있다.

paretodens.c의 코드

paretodens.c 파일 안에 우리가 paretodens 함수를 정의하였다고 가정하자. 그렇다면 paretodens.c는 다음의 코드를 포함하고 있어야 한다.

paretodens.c

#include <R.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* declare function */
void paretodens(double *x, double *alpha,
                double *x_m, int *mlength,
                double *density){
  
  // Intialization
  int n, i;
  n = *mlength;
  
  // For loop in C
  for (i = 0; i < n; i++){
    if (alpha[i] <= 0.0 || x_m[i] <= 0.0) {
      density[i] = NAN;
    } else {
      if (x[i] >= x_m[i]) {
        density[i] = log(alpha[i]) + 
          alpha[i] * log(x_m[i]) -
          (alpha[i] + 1.0) * log(x[i]);
      } else {
        density[i] = log(0.0);
      }
    }
  }
}

위의 코드를 보면 일반적인 C코드와 달리 앞에서 말한 #include <R.h>를 포함하고 있는 것을 확인할 수 있다. 자 이제, paretodens.c 파일을 컴파일한다. 여기서 기존의 C파일 컴파일과 다른 점은 C언어의 컴파일러를 사용하는 것이 아닌 R CMD SHLIB 명령어를 사용하여 함수를 라이브러리화 한다는 점이다. paretodens.c 파일이 있는 폴더를 기본 설정 폴더로 변경한 뒤, Rstudio 콘솔창에 다음의 명령어를 입력하자. 혹은 시스템 명령어 실행창에 따옴표안의 명령어를 입력한다.

system("R CMD SHLIB paretodens.c")

위의 명령어를 실행시키면, paretodens.c 파일있는 폴더에 paretodens.o파일과 paretodens.so 파일이 생겨난다. 우리는 paretodens.so 파일만 사용하게 된다(윈도우에서는 dll파일이 생성될 것이다). 이 .so파일은 shared object의 줄임말로 앞에서 C코드로 정의한 paretodens 함수의 정보를 담고있다고 생각하면 된다. 우리가 정의한 dpareto함수를 실행시키기 위해서는 paretodens 함수가 포함되어있는 .so 파일을 R로 불러와야 한다.

dyn.load("paretodens.so")

위의 명령어는 paretodens.so 파일의 내용을 R의 메모리 상에 올려주는 역할을 한다. 이 작업을 마친 후, 우리가 설정한 dpareto 함수를 실행시켜 보도록 하자.

all.equal(dpareto(3,1,2), 0.2222222222)
## [1] TRUE
all.equal(dpareto(1,3,2), 0.0)
## [1] TRUE
all.equal(dpareto(1,-3,2), NaN)
## Warning in dpareto(1, -3, 2): All parameters should be positive.
##  NaN is generated.
## [1] TRUE
all.equal(dpareto(1, 3,-2), NaN)
## Warning in dpareto(1, 3, -2): All parameters should be positive.
##  NaN is generated.
## [1] TRUE
all.equal(dpareto(3:5,1, 2), c(0.2222222222, 0.1250000, 0.0800000))
## [1] TRUE
all.equal(dpareto(6, 1, 2:4), c(0.05555555556, 0.08333333333, 0.11111111111))
## [1] TRUE
dpareto(3, 1, 2, log = TRUE) # -1.504077
## [1] -1.504077

이번에는 저번 [포스팅]에서 가장 좋은 성능을 보였던 Kyung Mo Kweon님의 dpareto4와 필자가 짠 dpareto3와의 속도를 비교해보도록 하자.

  • dpareto3 코드
dpareto3 <- function(x, alpha, x_m, log = FALSE) {
  mlength <- max(length(x), length(alpha), length(x_m))

  density <- vector(mode="numeric", length = mlength)
  x     <- rep(x, length.out = mlength)
  alpha <- rep(alpha, length.out = mlength)
  x_m   <- rep(x_m, length.out = mlength)

  #Filterling error condition & support condition
  invalid <- alpha <= 0 | x_m <= 0
  sub <- x >= x_m & alpha > 0 & x_m > 0

  # Invalid parameter
  if (sum(invalid) > 0){
    warning("All parameters should be positive.\n NaN is generated.\n")
  }
  
  #Give 'NaN' values for the inputs which contains errors.
  density[sub] <- log(alpha[sub]) + alpha[sub] * log(x_m[sub]) -
                  (alpha[sub] + 1) * log(x[sub])
  density[invalid] <- NaN
  density[!invalid & !sub] <- -Inf

  if (log) density
  else exp(density)
}
  • dpareto4 코드
dpareto4 <- function(x, alpha, x_m, log = FALSE) {
  density <- log(alpha) + alpha * log(x_m) - (alpha + 1) * log(x)
  density[x < x_m] <- -Inf
  density[alpha <= 0 | x_m <= 0] <- NaN
    
  if (any(is.nan(density))) {
    warning("All parameters should be positive.\n NaN is generated.\n")
  }

  if (log) density
  else exp(density)
}

성능비교 1차

library(microbenchmark)
set.seed(2017)
microbenchmark(
   dpareto(seq(1, 10, by = 0.01), 1:100, 3:5),
  dpareto3(seq(1, 10, by = 0.01), 1:100, 3:5),
  dpareto4(seq(1, 10, by = 0.01), 1:100, 3:5)
)
## Unit: microseconds
##                                         expr     min       lq     mean
##   dpareto(seq(1, 10, by = 0.01), 1:100, 3:5) 125.702 129.3765 167.0368
##  dpareto3(seq(1, 10, by = 0.01), 1:100, 3:5) 159.561 163.2875 331.3177
##  dpareto4(seq(1, 10, by = 0.01), 1:100, 3:5) 234.653 243.1230 268.8175
##    median      uq       max neval
##  132.8455 137.060  1865.149   100
##  165.8580 176.944 12140.573   100
##  249.4890 258.963  1636.426   100
all.equal(
  dpareto(seq(1,10, by = 0.01), 1:100, 3:5),
  dpareto3(seq(1,10, by = 0.01), 1:100, 3:5),
  dpareto4(seq(1,10, by = 0.01), 1:100, 3:5)
)
## [1] TRUE

C와 연결한 dpareto의 성능이 잘 나오는 것을 확인할 수 있다. 또한 한가지 의외의 결과가 나왔는데, 저번 포스팅에서 필자의 dpareto3보다 성능이 더 잘나왔던 dpareto4의 성능이 가장 느리게 나왔다. 따라서 필자는 다음의 비교를 한번 더 해보았다.

성능 비교 2차

set.seed(2017)
microbenchmark(
   dpareto(seq(1, 10, by = 0.001), 1:100, 3:5),
  dpareto3(seq(1, 10, by = 0.001), 1:100, 3:5),
  dpareto4(seq(1, 10, by = 0.001), 1:100, 3:5)
)
## Unit: microseconds
##                                          expr      min        lq      mean
##   dpareto(seq(1, 10, by = 0.001), 1:100, 3:5)  999.490 1021.7565 1677.5181
##  dpareto3(seq(1, 10, by = 0.001), 1:100, 3:5) 1272.560 1334.2535 1534.0877
##  dpareto4(seq(1, 10, by = 0.001), 1:100, 3:5)  819.375  851.6805  984.7511
##    median        uq       max neval
##  1043.937 1080.3965 47286.800   100
##  1354.818 1400.8750  2746.527   100
##   881.783  925.4835  2271.528   100
all.equal(
   dpareto(seq(1, 10, by = 0.001), 1:100, 3:5),
  dpareto3(seq(1, 10, by = 0.001), 1:100, 3:5),
  dpareto4(seq(1, 10, by = 0.001), 1:100, 3:5)
)
## [1] TRUE
dyn.unload("paretodens.so")

결과를 보면 이번에는 반대로 dpareto4의 성능이 오히려 C와 연결된 dpareto 보다 좋게 나온다.

개인적 의견

이러한 현상의 가장 큰 이유는 dpareto4의 경우는 recycling을 이용하여 필자가 선택한 입력값 만드는 방식을 피해서 짜여진 것에 있다고 생각한다. 하지만 이 방식의 경우 모든 함수값을 계산한 후 다시 계산한 케이스 중 0인 값과 NaN인 값을 중복해서 계산한다. 이러한 측면에서 필자가 짠 dpareto3의 경우는 케이스를 나눈 후 한번씩만 계산하는 방식이어서 효율성이 있다. 이러한 이유로 dpareto3가 입력값을 만드는 시간이 비교적 짧은 1차 비교의 경우 dpareto4 보다 더 효율적이다.

하지만 dpareto4는 2차 비교처럼 입력값의 벡터가 커지면 커질 수 록 dpareto4의 recycling의 효율성이 증가해서 중복 계산의 비효율성을 커버할 수 있다는 게 장점이다. 필자가 작성한 dpareto의 경우도 dpareto3와 똑같은 입력 벡터를 만드는 과정을 취하고 있는데, 이러한 점에서 dpareto3와 똑같은 느려짐 현상이 발생한다고 생각한다. dpareto가 dpareto4 함수보다 모든경우의 성능 테스트에서 우월해지기 위해서는 입력벡터 만드는 방식 전체를 C안에서 해버리는게 좋다. 하지만 이것은 사용자의 가치판단의 문제가 생기는데, 입력벡터를 만드는 작업까지 C에서 해버린다면 그건 그냥 C 프로그래밍이지 R프로그래밍은 아니라는게 필자의 견해이다. (하지만 만약 R 패키지 제작을 목적으로 한다면 이야기가 살짝 달라지긴한다.)

  • 요약
  1. 최적화된 R코드의 경우, 필자가 생각했던 것보다 빠름.
  2. 정말 빠르게 하고싶다면 하이브리드도 생각해보면 좋을것 같음.
  3. 필자의 코드와 dpareto4의 장점을 합친 보다 최적의 코드가 있을것 같기도 함.


SHARE TO

신고


티스토리 툴바