隙間

こちらはいわゆる仕事用のもの.

Call C in your R

この記事は明治大学 Advent Calendar 2018 - Qiitaの13日目の記事です.

修論の合間を縫ってめちゃくちゃ急いで書きました.不備があったら,ぜひ教えてください.可能なかぎり修正します.


最近,僕のインターネットの知り合いの方の中で洋画を夜な夜な漁るのが流行っていて, その界隈の中で「君の名で僕を呼んで」 という映画がとにかく良い,という話を耳にするのです.その英題が ”call me by your name.”

なので,今回はそれにちなんで「Rスクリプトの中からCの実行ファイルを実行し,計算速度を上げるtips」としてpipeと呼ばれる技術を紹介します.

現在では後発でRcppという大変便利なライブラリが開発され,ここで紹介するより高速に処理をすることができる技術もありますが,呼び出すのがC++で,タイトルの語感を崩したくはないので,  pipe自体がR&Cという言語の組み合わせに限らず,多くの言語で実装可能な汎用性の高い技術なので,こちらを紹介します.

Call C in your R.



pipeとは

pipeとは,ある実行ファイルが処理(process)を進めている最中に,途中まで処理をした(or 自身が入力として受け取ったデータ)を別の実行ファイルに入力として渡し,その実行ファイル内で設計された演算処理をしてもらったあとにデータを出力として受け取る作業のことを言います.

CやPythonに代表される手続き型言語において,プログラミングの基本概念の1つである関数でも,関数はメインモジュール(Cでのmain関数に相当)から引数(argument)を入力として受け取り,演算処理を施した後,返却値(return value)として出力するというのが基本の流れです.また,C言語でのscanf()printf()関数に代表される標準入出力(standard input and output)もユーザー(人間)とコンピューター間の入力と出力を担っているものですよね.

このようにコンピューターにまつわる情報処理では様々な階層で「入力→演算→出力」が行われています.今回紹介するpipeもまた,「実行ファイルと実行ファイル間の入力→演算→出力」を実装する手段だと考えてもらうと,後述の内容が理解しやすいでしょう.

実行時間・処理速度・計算量

実行時間

また今回の記事はプログラムの速さ・処理能力といった項目を扱うので,簡単にですが,これらを議論・検証するための用語を整理しておきます.概念的には極めて重要なものですが,それと同時に直感的に理解しづらい部分もあるため,1度最後まで読んだ後に時間をおいて,ここに帰ってきてもらっても大丈夫です.

これ以降,ある処理を完了させるまでの時間を実行時間 (running time)と定義します.教科書や分野によって,用語の用法に差がありますが,今回はこのように定義をすることにします.後述する「処理速度」「計算量」についても同様です.*1

処理速度

処理速度とは,実行ファイルがプログラム内に記述された入力・演算・出力等の処理を完了させるまでの速さを指します.

一般論として,同一の環境で同一の処理を行うとき,同一のアルゴリズムを実装しても,実装に用いたプログラミング言語によってその処理の実行時間には差があります.これをプログラミングの世界では「処理速度が違う」と表現することが多いです.メモリへのアクセスの仕組みや,メモリでのデータ保持の方法の違いなど処理速度の違いの原因となるものは様々あります.よく「XXは速いよね」とか「XXXはfor文遅いんだよねぇ」と言われるのは往々にしてこの処理速度のことを指しています.

計算量

前述した処理速度の他に,実行時間を決める要因の一つに計算量と呼ばれる概念が存在します. 計算量(computational complexity)とは計算可能な問題において,その計算を解くためのを完了させるまでの計算資源の必要量を指します.

より厳密には,計算を完了させるための演算(四則演算・真偽判定などなど,プログラミングにおけるステップ)回数を評価する時間計算量とその計算を完了させるまでに必要な情報の記憶領域の規模を表す空間計算量に分けて考えます.*2

時間計算量は,ランダウの記号O(\cdot)と処理の対象となるデータの入力サイズnを用いて,n多項式時間として表現します.アルゴリズムの良し悪しを評価する際の重要な指標の一つです. 詳しくは以下のページを参照してください:


qiita.com


時間計算量の概念を知っていると,いくら,処理速度が早くても,処理にかかる時間計算量が大きければ,結果的に実行時間が遅くなるというのが容易に想像がつくでしょう.このように実行時間に関わるのは時間計算量です.以降,断りがない場合は計算量とは時間計算量のことを指すこととします.

ベンチマーク

プログラムにおいて,実際に実装した処理の実行時間を定量的に評価するために,今回は検証する2言語(RとC)を用いて,同一の処理を同一のアルゴリズムで実装し,処理速度を評価します.

今回検証で使うのはこの手のテストでおなじみのフィボナッチ数列の値を求める計算です:

F_{n+2} = F_{n+1} + F_{n}(n \geqq 0),

F_{1} = 1, F_{0} = 1.

また,検証に用いた環境は以下の通りです:

macOS 10.12.6
MacBook Pro (Retina, 13-inch, Early 2015)
2.7 GHz Intel Core i5
8 GB

また,C言語,Rの各スクリプト内で実行時間を測定するために,それぞれ以下のライブラリ・モジュールを使用しました. それぞれの記事を参考にしてください:


stackoverflow.com


hiratake55.hatenadiary.jp


さらに,今回の検証ではフィボナッチ数列の計算を実装する際にメモ化を用いた実装をしています. メモ化について知らない人は後述のメモ化の項目をあらかじめ読んでください.次以降の項目で,実際にベンチマークの結果を見ていきましょう.

(補足)メモ化

フィボナッチ数列は数学的な定義に基づけば,再帰的(recursively)に計算することができます.具体的にどのようにコンピューター内部で計算されているかはこのアニメーション などを参考にしてください.

しかし,プログラムで再帰計算をそのまま実装してしまうと,計算時間が指数関数的に増加することが一般的に知られています.*3 *4

おそらく,前述したアニメーションを見れば,この処理速度の低下(実行時間の増加)は項数が増えるほど,再起の深さが大きくなることに起因することが直感的に理解できるでしょう.

そこで,このような再帰処理を実装する際にプログラミングでは(とりわけ,競技プログラミングの分野)ではメモ化と呼ばれる技術を用いることで,再帰深度を抑え,計算量を抑えることがしばしば行われます.

詳しくは以下のページを参照してください:


www.slideshare.net


qiita.com


Fibonacci in C

まずはC言語単独でフィボナッチ数列の計算のベンチマークを行います.実際のコードは以下の通りです.

implementation in C
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

#define N /* input num of terms in the sequence here. */

int Fibonacci(int i);
int Fmemo[N];

int main(int argc, const char * argv[]) {
    int long a=0;
    struct timeval tval_before, tval_after, tval_result;
    
    int i;
    for (i = 0; i < N; i++) {
        Fmemo[i] = 0;
    }
    
    Fmemo[0] = 1;
    Fmemo[1] = 1;
   
    gettimeofday(&tval_before, NULL);
    a = Fibonacci(N);
    gettimeofday(&tval_after, NULL);
    timersub(&tval_after, &tval_before, &tval_result);
    printf("Time elapsed: %ld.%09ld\n", (long int)tval_result.tv_sec, (long int)tval_result.tv_usec);
       
    return 0;
}

int Fibonacci(int i){
    switch (i) {
        case 0:
            return 1;
            break;
            
        case 1:
            return 1;
            break;
            
        default:
            
            if (Fmemo[i]!=0) {
                return Fmemo[i];
            }else{
                Fmemo[i] = Fibonacci(i-1) + Fibonacci(i-2);
                return Fmemo[i];
            }
     
            break;
    }
}
result

結果は以下の通り*5

N time (sec)
20 1.0E-10
200 4.0E-9
2000 4.6E-8
20000 8.14E-7

Fibonacci in R

次にR単独でフィボナッチ数列の計算のベンチマークを行います.実際のコードは以下の通りです.

implementation in R
N <- ## input num of terms in the sequence here. ##
Fmemo <- rep(0,N)
Fmemo[1] <- 1
Fmemo[2] <- 1

Fibonacci <- function(i){
  if(i == 1 || i == 2){
    return(1)
  }else{
    if(Fmemo[i]!=0){
      return(Fmemo[i])
    }else{
      Fmemo[i] <- Fibonacci(i-1) + Fibonacci(i-2)
      return(Fmemo[i])
    }
  }
}

system.time(
  print(Fibonacci(N))
  )
result

結果は以下の通り:

N time (sec)
20 3.0E-2
200 356+
2000 ---
20000 ---

Rの場合に関しては,200\leqq Nの場合は計算が終了する予感もしないため検証を断念しました... (760secくらいまでは待った.)

Call C in your R

ここまでで検証したように,Rでそのまま数列の計算を実装すると,実行時間をとても多く要します.そもそも,このような用途にRが適していないため,言語の設計構造上不向きであることや,ルーチン(関数)の呼び出しに時間がかかることが原因のようです.

それではこの問題を解決するために,Nの値を入力し,フィボナッチ数列の値を出力する関数を.cファイルで記述し,その関数をRのネイティブ関数である.Cインターフェースを介して,サブルーチンとして受け取る方法を実装します.

anythingbutrbitrary.blogspot.com

具体的な手順は以下の通りです:

  1. C言語で関数を実装

  2. .cファイルを$ R CMD SHLIBコンパイル

  3. Rスクリプトを実行

順を追って確認します.


C言語で関数を実装

Rの .Cインターフェースを介して入力を受け取るとき,.cファイル内で実装された関数は引数を ポインタ変数として受け取ります.したがって,関数を実装する際はポインタ型変数を受け取る関数として実装します.また,R側で呼び出す関数は返却値を持たないvoid 型でなければいけません.

この仕様の関係上,今回は項数(N)の値を持つ変数の上にF_Nの値を上書きするような方法で実装されていることに注意してください.

さらに今回の場合,メモ化に使う配列を用いる関係で以下のように実装します.(staticC とか使えばもうちょい賢くできそう...)

fibo_routine.c

long int Fmemo[30000];
long int Fib(int i);

void Fibonacci(int *N){
    int i;
    for(i = 0; i < 30000; i++){
        Fmemo[i] = 0;
    }
    
    Fmemo[0] = 1;
    Fmemo[1] = 1;
    
    *N = Fib(*N);
}

long int Fib(int i){
    switch (i-1) {
        case 0:
            return 1;
            break;
            
        case 1:
            return 1;
            break;
            
        default:
            
            if (Fmemo[i]!=0) {
                return Fmemo[i];
            }else{
                Fmemo[i] = Fib(i-1) + Fib(i-2);
                return Fmemo[i];
            }
            
            break;
    }
}
.cファイルを$ R CMD SHLIBコンパイル

次に fibo_routine.cコンパイルしてバイナリ(実行ファイル)の形にします.ここで通常であれば,cc コマンド等でコンパイルしますが,今回は以下のように:

$ R CMD SHLIB fibo_routine.c

コンパイルします.

正常にコンパイルできれば,fibo_routine.o fibo_routine.so というバイナリファイルが生成されるはずです.

Rスクリプトを実行

前段階で生成されたバイナリファイルのうち fibo_routine.so 側を使用し,R側で呼び込みます.

dyn.load("fibo_routine.so") の1行でバイナリファイルを読み込み,次の.C()の部分で.cファイル内で実装したFibonacci 関数を呼び出しています.

dyn.load("fibo_routine.so")
.C("Fibonacci", N = as.integer(##input the val of N here##))

実際にスクリプトを実行した結果は以下の通りです:

> dyn.load("fibo_routine.so")
> .C("Fibonacci", N = as.integer(20))
$N
[1] 6765

 F_{20} = 6765が出力されていることを確認してください.

比較結果

それでは以上の結果を踏まえて,R単独での実装と.C込みでのベンチマークを比較します.

小数点以下を正確に測定するためにベンチマークの方法を若干変更しています.*6

dyn.load("fibo_routine.so")

options(digits.secs = 6) # This is set so that milliseconds are displayed
start.time <- Sys.time()
system.time(.C("Fibonacci", N = as.integer(20000)))
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken

結果は以下の通り:

N R standalone (sec) R w/ .C interface (sec)
20 3.0E-2 1.04E-1
200 356+ 1.12E-1
2000 --- 1.14E-1
20000 --- 1.25E-1

最初の N=20の時こそ,効果は見られませんが,それ以降を比較すると差は歴然です.

早いね.

こうすることで,シミュレーション本体はCで実行し,結果の解析をRでやるというフローがワンフローでできるなど恩恵は様々です.

興味があれば,冒頭で紹介したRcppも活用して見てください.そして,残り2枠のアドベントカレンダーをぜひ埋めましょう.

(完)

*1:ここらへんの計算複雑性の理論は普通に怪しいので,もし,間違ってたらコメントください.

*2:チューリングマシンでは時間計算量が「動作ステップ数」空間計算量が「テープ長」に相当.詳しくは現象数理学科「情報処理」「数学の方法」「アルゴリズム論」他

*3:http://www.cs.tsukuba.ac.jp/~kam/lecture/gairon1/SS1-2013-algorithm.pdf

*4:厳密には,  O\left( \frac{1+\sqrt{5}}{2}\right)^n

*5:https://www.google.co.jp/search?q=1.0e-10%E3%81%A8%E3%81%AF&rlz=1C5CHFA_enJP709JP709&oq=1.0E-10&aqs=chrome.1.69i57j0l5.5343j0j4&sourceid=chrome&ie=UTF-8

*6:https://stackoverflow.com/questions/7546946/get-execution-time-in-milliseconds-in-r/36646392