본문으로 바로가기



Make your own R package! - 8. C언어에서 recycling 구현과 functional programming

category R Programming/Rpackage 2017.08.05 19:00

이전 포스팅에서 우리는 R함수와 C언어 짜여진 함수를 연결하는 방법에 대하여 알아보았다. 이전 포스팅에서는 R 함수를 이용하여 입력값의 길이가 다른 경우 lengthout 옵션을 이용하여 같은 길이의 벡터로 바꿔주었는데, 이 부분을 이번에는 C에서 구현해보도록 하자.

C 함수에서의 recyling 구현

C 함수에서 R에서의 recycling을 구현하기 위해서는 % 연산자를 잘 사용하면 된다. 필자가 정의한 함수의 코드를 살펴보도록 하자. 아래의 코드에서 입력값들은 x, alpha, x_m 벡터와 각각의 길이를 담고있는 input_len 벡터, 그리고 가장 긴 벡터의 길이 max_len, log를 처리하기 위한 indicator, 결과를 전달받을 density 벡터를 의미한다.

가장 핵심이 되는 부분은 for 구문인데 i는 0에서 부터 n-1까지 움직이며(n은 세 개 입력 벡터 중 가장 긴 벡터의 길이), 각 i에 대응하는 값을 [i % length of each vector]로 찾아서 대입하는 방식을 사용한다. 잘 이해가 되지않는 독자들은 종이와 펜을 들고 직접 써가면서 이해하면 훨씬 쉬울것이라 생각한다.

void paretodens1(double *x, double *alpha, double *x_m, 
                int *input_len, int *max_len, int *is_log,
                double *density){
  
  // Intialization
  int i, n, x_len, alpha_len, x_m_len, IsLog;

  // Prepare recycling
  n = *max_len;
  IsLog = *is_log;
  x_len = input_len[0];
  alpha_len = input_len[1];
  x_m_len = input_len[2];
  
  // For loop in C
  for (i = 0; i < n; i++){
    if (alpha[i % alpha_len] <= 0.0 || x_m[i % x_m_len] <= 0.0) {
      density[i] = R_NaN;
    } else {
      if (x[i % x_len] >= x_m[i % x_m_len]) {
        density[i] = log(alpha[i % alpha_len]) + 
            alpha[i % alpha_len] * log(x_m[i % x_m_len]) -
                (alpha[i % alpha_len] + 1.0) * log(x[i % x_len]);
      } else {
        density[i] = R_NegInf;
      }
    }
  }
  
  if (!IsLog) {
      for ( i = 0; i < n; i++){
          density[i] = exp( density[i] );
      }
  }
}

필자의 위의 코드가 약간 가독성이 떨어지는 것 같아서, for 루프를 들어가기전 벡터를 미리 만드는 방법을 택하는 형식의 코드를 하나 더 정의하였다. 아래의 형식을 택하면 코드의 길이가 길어지고, 입력벡터를 하나 더 만드는 셈이라서 비효율적이지만, 코드의 가독성은 증가한다.

void paretodens5(double *x, double *alpha, double *x_m, 
                 int *input_len, int *max_len, int *is_log,
                 double *density){
  
  // Intialization
  int i, n, x_len, alpha_len, x_m_len, IsLog;
  
  // Prepare recycling
  n = *max_len;
  IsLog = *is_log;
  x_len = input_len[0];
  alpha_len = input_len[1];
  x_m_len = input_len[2];
  
  double my_x[n], my_alpha[n], my_x_m[n];
  
  // Recyling
  if (x_len == n){
    for (i = 0; i < n; i++){
      my_x[i] = x[i];
      my_alpha[i] = alpha[i % alpha_len];
      my_x_m[i] = x_m[i % x_m_len];
    }
  } else if (alpha_len == n) {
    for (i = 0; i < n; i++){
      my_x[i] = x[i % x_len];
      my_alpha[i] = alpha[i];
      my_x_m[i] = x_m[i % x_m_len];
    }
  } else {
    for (i = 0; i < n; i++){
      my_x[i] = x[i % x_len];
      my_alpha[i] = alpha[i % alpha_len];
      my_x_m[i] = x_m[i];
    }
  }
  
  // For loop in C
  for (i = 0; i < n; i++){
    density[i] = R_NegInf;
  }
  
  for (i = 0; i < n; i++){
    if (my_alpha[i] <= 0.0 || my_x_m[i] <= 0.0) {
      density[i] = R_NaN;
    } else if(my_x[i] >= my_x_m[i]){
      density[i] = log(my_alpha[i]) + 
        my_alpha[i] * log(my_x_m[i]) -
        (my_alpha[i] + 1.0) * log(my_x[i]);
    }
  }
  
  // Log
  if (!IsLog) {
    for ( i = 0; i < n; i++){
      density[i] = exp( density[i] );
    }
  }
}

다음은 필자가 정의해놓은 C코드를 담고 있는 파일이다. 아래의 파일을 이용하여 so파일 혹은 dll파일을 만들도록 하자.

paretodens.c

system("R CMD SHLIB paretodens.c")
  • 한 가지 C코드를 짤 때 주의할 사항은 R과 연동하여 정의하는 C코드의 경우 #include <R.h>를 넣어주고, NA_REAL, R_NaN, R_PosInf or R_NegInf 등을 사용할 것을 권장하고 있다는 점이다.

Functional programming을 사용한 .C 파일 정의하기

위의 paretodens.c 파일에는 6개의 함수가 들어있다. recyling을 구현하기 위해 필자가 여러개의 C코드를 짜면서 깨달은 것은 필자의 C함수들이 R로 부터 항상 같은 형식의 입력값을 받는다는 것이었다. 즉, C로 넘겨주는 입력값은 똑같은데, 값을 받고난 후의 C 함수는 다른 상황이어서 이 부분을 좀더 효율적으로 짜보기로 하였다. 따라서 .C 파일을 사용한 R함수는 다음과 같은 함수를 만들어주는 함수로 정의하여 놓았다. 이렇게 정의해놓으면 적용하는 C함수의 이름만 바꿔주면 각기 다른 함수를 쉽게 정의할 수 있다.

make_dpareto <- function(fname = "paretodens") {
    function(x, alpha, x_m, log = FALSE){
        input_len <- c(length(x), length(alpha), length(x_m))
        max_len <- max(input_len)
        
        # Input vectorizing
        density <- vector(mode="numeric", length = max_len)
        
        # Call .C
        result <- .C(fname,
                     as.double(x),
                     as.double(alpha),
                     as.double(x_m),
                     as.integer(input_len),
                     as.integer(max_len),
                     as.integer(log),
                     value = as.double(density))$value
        if (any(is.nan(result))) {
          warning("\n All parameters should be positive.\n NaN is generated.")
        }
        result
    }
}

위의 make_paretodens 함수의 경우, 입력값과 같은 이름의 C함수를 찾아서 적용시켜주는 함수를 만들어주게 된다. 즉, 우리가 paretodens1, paretodense2, 등과 같은 이름의 함수를 정의해놓았다면, 다음과 같이 각기 다른 함수를 효율적으로 정의할 수 있도록 해주는 것이다.

# These dpareto1~5 pass everything to C
dpareto1 <- make_dpareto("paretodens1")
dpareto2 <- make_dpareto("paretodens2")
dpareto3 <- make_dpareto("paretodens3")
dpareto4 <- make_dpareto("paretodens4")
dpareto5 <- make_dpareto("paretodens5")

함수 테스트

먼저 순수 R 함수로 짠 dpareto_pure를 다음과 같이 정의한다.

# dpareto_pure is pure R code dpareto
dpareto_pure <- function(x, alpha, x_m, log = FALSE) {
  density <- suppressWarnings(log(alpha) + alpha * log(x_m) - (alpha + 1) * log(x))
  density[suppressWarnings(x < x_m)] <- -Inf
  density[suppressWarnings(alpha <= 0 | x_m <= 0)] <- NaN
  
  if (any(is.nan(density))) {
    warning("\n All parameters should be positive.\n NaN is generated.")
  }
  
  if (log) density
  else exp(density)
}

위의 함수와 우리가 짠 함수가 같은 값을 가지는가에 대하여 다음과 같이 테스트해보자.

dyn.load("paretodens.so")
    all.equal(dpareto1(1:4, 2:4, 3:8),
              dpareto_pure(1:4, 2:4, 3:8))
## [1] TRUE
    all.equal(dpareto2(3, -2:4, 3:5),
              dpareto_pure(3, -2:4, 3:5))
## Warning in dpareto2(3, -2:4, 3:5): 
##  All parameters should be positive.
##  NaN is generated.
## Warning in dpareto_pure(3, -2:4, 3:5): 
##  All parameters should be positive.
##  NaN is generated.
## [1] TRUE
    all.equal(dpareto3(4, -2:4, 1:7),
              dpareto_pure(4, -2:4, 1:7))
## Warning in dpareto3(4, -2:4, 1:7): 
##  All parameters should be positive.
##  NaN is generated.
## Warning in dpareto_pure(4, -2:4, 1:7): 
##  All parameters should be positive.
##  NaN is generated.
## [1] TRUE
    all.equal(dpareto1(3,1,2), 0.2222222222)
## [1] TRUE
    all.equal(dpareto1(1,3,2), 0.0)
## [1] TRUE
    all.equal(dpareto2(1,-3,2), NaN)
## Warning in dpareto2(1, -3, 2): 
##  All parameters should be positive.
##  NaN is generated.
## [1] TRUE
    all.equal(dpareto2(1, 3,-2), NaN)
## Warning in dpareto2(1, 3, -2): 
##  All parameters should be positive.
##  NaN is generated.
## [1] TRUE
    all.equal(dpareto3(3:5,1, 2),
              c(0.2222222222, 0.1250000, 0.0800000))
## [1] TRUE
    all.equal(dpareto3(6, 1, 2:4),
              c(0.05555555556, 0.08333333333, 0.11111111111))
## [1] TRUE
    dpareto4(3, 1, 2, log = TRUE) # -1.504077
## [1] -1.504077
dyn.unload("paretodens.so")

성능 테스트

짠 함수의 성능을 비교해 보기 위하여 다음과 같은 코드로 테스트 하여보았다.

dyn.load("paretodens.so")
    library(microbenchmark)
    test_performance <- function(n, seed = 2017){
      set.seed(seed)
      my_seq <- seq(1, 10, by = 10^(-n))
      microbenchmark(
          dpareto_pure(my_seq, 1:10, 3:7),
              dpareto1(my_seq, 1:10, 3:7),
              dpareto2(my_seq, 1:10, 3:7),
              dpareto3(my_seq, 1:10, 3:7),
              dpareto4(my_seq, 1:10, 3:7),
              dpareto5(my_seq, 1:10, 3:7)
      )
    }
    test_performance(1)
## Unit: microseconds
##                             expr     min       lq      mean   median
##  dpareto_pure(my_seq, 1:10, 3:7) 137.482 141.9020 146.38680 144.8990
##      dpareto1(my_seq, 1:10, 3:7)  15.323  16.1710  17.34816  16.8105
##      dpareto2(my_seq, 1:10, 3:7)  15.203  16.2905  17.41094  16.7205
##      dpareto3(my_seq, 1:10, 3:7)  18.913  20.3545  21.36990  20.8920
##      dpareto4(my_seq, 1:10, 3:7)  17.736  18.8160  19.92723  19.4310
##      dpareto5(my_seq, 1:10, 3:7)  15.215  16.1315  16.89356  16.4715
##        uq     max neval
##  148.2690 252.771   100
##   18.3870  22.012   100
##   17.2945  50.764   100
##   21.8395  27.297   100
##   20.4190  24.691   100
##   17.3010  20.990   100
    test_performance(2)
## Unit: microseconds
##                             expr     min       lq      mean   median
##  dpareto_pure(my_seq, 1:10, 3:7) 188.172 194.4670 200.83990 198.7280
##      dpareto1(my_seq, 1:10, 3:7)  78.403  79.5225  81.85638  81.3600
##      dpareto2(my_seq, 1:10, 3:7)  79.015  79.8960  82.26114  81.0895
##      dpareto3(my_seq, 1:10, 3:7) 117.547 119.1865 121.99164 120.8690
##      dpareto4(my_seq, 1:10, 3:7) 104.807 106.3275 125.71228 107.8960
##      dpareto5(my_seq, 1:10, 3:7)  76.513  77.5120  79.40653  78.5090
##        uq      max neval
##  204.5315  301.891   100
##   83.5210   90.191   100
##   83.6980  119.265   100
##  123.9725  133.501   100
##  110.8440 1809.466   100
##   81.0235   95.208   100
    test_performance(3)
## Unit: microseconds
##                             expr     min        lq      mean    median
##  dpareto_pure(my_seq, 1:10, 3:7) 617.803  689.4115  758.4695  703.1955
##      dpareto1(my_seq, 1:10, 3:7) 625.965  689.9620  721.8309  706.2040
##      dpareto2(my_seq, 1:10, 3:7) 643.389  707.6945  715.3622  710.6275
##      dpareto3(my_seq, 1:10, 3:7) 987.133 1092.3020 1115.4535 1098.5250
##      dpareto4(my_seq, 1:10, 3:7) 904.662  975.9680 1047.0937  980.4440
##      dpareto5(my_seq, 1:10, 3:7) 624.566  683.1225  719.5773  688.3575
##         uq      max neval
##   748.5895 2263.076   100
##   713.5195 2175.739   100
##   716.0030  798.512   100
##  1107.4885 2459.230   100
##   991.3140 2629.454   100
##   695.2530 2193.391   100
    test_performance(4)
## Unit: milliseconds
##                             expr      min       lq      mean   median
##  dpareto_pure(my_seq, 1:10, 3:7) 4.586413 4.736482  5.447596 5.299397
##      dpareto1(my_seq, 1:10, 3:7) 5.695688 5.821592  7.233784 6.597784
##      dpareto2(my_seq, 1:10, 3:7) 5.721037 5.842486  6.443545 6.129651
##      dpareto3(my_seq, 1:10, 3:7) 8.888376 9.030519  9.804542 9.176205
##      dpareto4(my_seq, 1:10, 3:7) 7.946121 8.095397 10.173255 8.562967
##      dpareto5(my_seq, 1:10, 3:7) 5.555754 5.662800  8.264295 5.925059
##         uq       max neval
##   6.040546  7.211671   100
##   7.126734 71.597985   100
##   7.084509  8.351595   100
##  10.464908 12.342867   100
##   9.475160 74.645336   100
##   6.874781 72.521525   100
dyn.unload("paretodens.so")

흥미로운 사실은 입력 벡터값의 길이가 늘어나면 늘어날수록 순수 R코드 함수의 성능과 .C 방법을 사용한 함수들간의 차이가 좁혀진다는 것이다. 필자의 개인적인 견해는 다음과 같다. (어디까지나 추측이다.)

.C방법은 R과 C의 다리 역할을 해주는데, 이 과정에는 데이터를 R에서 C로 넘겨주는 과정이 포함되어 있다. 따라서 R에서 C로 넘겨주는 벡터 크기가 커지면 커질수록 두 프로그램간 데이터 전송시간도 함수의 실행 시간에 추가되어 지는 것이다. 우리가 짠 함수의 경우 입력 벡터의 길이가 10이면 계산을 마친 후 결과 벡터도 10인 벡터가 된다. 따라서 데이터의 전송시간과 결과값을 받는 시간이 똑같다. 이렇게 C로 전송하는 벡터 대비 전송받는 벡터의 크기가 같거나 큰 경우는 두 프로그램 간 데이터 전송 시간이 프로세스의 병목구간이 되어버릴 가능성이 있다. 이제까지 말한 것은 사실 필자의 추측인데, 이것을 확인하기 위해서는 전속벡터 대비 결과 값이 상수 하나 같은, 예를 들어 평균이나 분산을 구하는 프로그램을 짜서 결과를 비교해보면 될 것이라 생각한다. 이 부분은 계속 포스팅을 진행해 나가면서 확인해보려 한다.

이 가정이 맞다면 결과적으로 C 혹은 다른 R외의 다른 언어의 힘을 빌릴지 말지는 자신이 어떤 함수를 짜느냐에 따라 달라질 수 있다. 이러한 과정을 알고 잘 함수를 설계한다면 훨씬 효율적인 프로그램이 탄생할게 될 것이라 생각한다.


SHARE TO

신고


티스토리 툴바