Rcppで定義するC++の関数に関数を渡すときはRcpp::FuncでOK

この辺

の例では1引数しか渡さないケースばかりだったので、複数引数の関数を渡せることに気がつかなかった。以下では適当な関数(someFunc)と、Rから関数を引数に取って処理を行う関数callFuncを定義している。

library(Rcpp)
sourceCpp(code='
  #include <Rcpp.h>
  using namespace Rcpp;

  // [[Rcpp::export]]
  NumericVector someFunc(NumericVector x, double t, CharacterVector c) {
    return (x+t);
  }
  // [[Rcpp::export]]
  NumericVector callFunc(NumericVector x, double t, CharacterVector c, Rcpp::Function func) {
    return func(x, t, c);
  }
')

実際に使ってみると・・・

> f <- function(x, t, s){x+t}
> callFunc(1:10, 1, "A", f)
 [1]  2  3  4  5  6  7  8  9 10 11
> callFunc(1:10, 1, "A", someFunc)
 [1]  2  3  4  5  6  7  8  9 10 11

となる。第三引数のCharacterVectorは、ちゃんと渡せるかを確認したかったただのおまけ。

参考

Rcpp::Funcの実装はここ

Rcpp使ってパフォーマンスを測ったらに分類される…?

こういうどうでもいい関数があったとする。

hoge <- function(x) {
  result = x + 1;
  for(i in 1:10000000){
    result <-  result + 1;
  }
  result
}

この関数の実行速度+パフォーマンスを測るにはRprof関数を使って以下のように書く。

> Rprof(tmp <- tempfile())
> a <- hoge(1:10)
> Rprof()
> summaryRprof(tmp)
$by.self
       self.time self.pct total.time total.pct
"hoge"      2.76    61.88       4.46    100.00
"+"         0.88    19.73       0.88     19.73
":"         0.82    18.39       0.82     18.39

$by.total
       total.time total.pct self.time self.pct
"hoge"       4.46    100.00      2.76    61.88
"+"          0.88     19.73      0.88    19.73
":"          0.82     18.39      0.82    18.39

$sample.interval
[1] 0.02

$sampling.time
[1] 4.46

んで、これが遅いなーと、計算のボトルネックだな―と思ってRcppを使って書きなおすと以下のようになる。

library(Rcpp)
sourceCpp(code='
  #include <Rcpp.h>
  using namespace Rcpp;
  
  // [[Rcpp::export]]
  NumericVector hoge(NumericVector x) {
    NumericVector result = x + 1;
    for(int i = 0; i < 10000000; i++){
      result = result + 1;
    }
    return result;
  }'
)

どうようにパフォーマンスを測ると

> Rprof(tmp <- tempfile())
> a <- hoge(1:10)
> Rprof()
> summaryRprof(tmp)
$by.self
              self.time self.pct total.time total.pct
"<Anonymous>"      0.14      100       0.14       100

$by.total
              total.time total.pct self.time self.pct
"<Anonymous>"       0.14       100      0.14      100
"hoge"              0.14       100      0.00        0

$sample.interval
[1] 0.02

$sampling.time
[1] 0.14

たしかに高速化されている。高速化されて速くなったのはいいんだけど、Anonymousって・・・C++を呼ぶ際の内部関数の名前なんだろうか。要調査。

追記

この辺と同じか。中で直に関数のポインタを呼んでるときは、Anonymous扱いされる模様。

Hadley、パフォーマンス計測用のパッケージ、複数作ってるっぽい。

最終コミットから察するにlineprofを伸ばしていくっぽい。そういえば昔試していた。

呼び出し関係を可視化するproftoolsなんてのもあるのか。

金融商品を記述するためのDSL(Domain Specific Language, ドメイン特化言語)っぽいもんをC++でやってみた(※元ネタ『Composing Contracts: An Adventure in Financial Engineering』)

はじめに


の6章「金融取引契約の書き方」をC++でオレオレ実装するとどうなるのか?って話をしましょうという話を某氏としたのでやってみんとす。C++辛い、F#とかがよかった。

これは一体何?

金融商品を記述するためのDSL(Domain Specific Language, ドメイン特化言語)を作って、それを実際評価(evaluation)して金融商品のお値段を計算してみましょうという枠組み?のお話。元の本とかだとHaskellなんかで実装するのがセオリーっぽい。自分の呟きでアレですが、Bloombergも既にOCamlで実装した模様。

あと、Githubの方にScala実装を上げてくれてた方がいた*1

オレオレ実装

割引債とヨーロピアンで力尽きました。C++辛い、lazyないの辛い、いちいちsmart pointer辛い、全体的に辛い。main部分は以下。それっぽくは見えるし、ここまでの答えは合ってるようだ。その他の実装について正しいかはわからん。わからんが、Githubにはアップしておいた。

int main()
{
	//One unit
	{
		ContractPtr c = one(JPY);
		eval_and_show(c);
	}
	//Zero Coupon bond
	{
		ContractPtr c = zcb(mkDate(365), 10.0, JPY);
		eval_and_show(c);
	}
	//Zero Coupon bond(port)
	{
		ContractPtr c1 = zcb(mkDate(365), 10.0, JPY);
		ContractPtr c2 = zcb(mkDate(365 * 2), 10.0, JPY);
		eval_and_show(and(c1, c2));
	}
	//European option
	{
		ContractPtr c = european(mkDate(365), zcb(mkDate(365 * 2), 10.0, JPY));
		eval_and_show(c);
	}
	return 0;
}

*1:当記事下部、参考を参照

カジュアルにC++11 in R

以下のプラグイン

plugins=c("cpp11")

をcppFunctionに指定するだけ。Winでもいけた。C++11なんで戻り値をdecltypeにしてやろうと思ったら、できなかった。windowsgccが古いせいかな?

library(Rcpp)
cppFunction('
  std::vector<double> twoTimes(std::vector<double> xs) 
  {
    for(auto &x : xs)
    {
      x *= 2.0;
    }
    return xs;
  }', plugins=c("cpp11"))

こんなかんじで。

> twoTimes(1:10)
 [1]  2  4  6  8 10 12 14 16 18 20

こうなる。

意外に速いRタソ compared to C++

ちょいとRcppをがっつり使ってみようと、その試行錯誤記録が続く予定。
速度検証のテストコードとして、ランダムウォークする系列データを生成するコードを書いた。

library(Rcpp)
sourceCpp(code='
  #include<vector>
  #include <functional>
  #include <numeric>
  #include <Rcpp.h>

  using namespace Rcpp;
  // [[Rcpp::export]]
  NumericVector randomWalk(int size)
  {
    RNGScope scope;
    NumericVector rand(Rcpp::rnorm(size));
    NumericVector cum(Rcpp::rnorm(size));
    std::partial_sum(rand.begin(), rand.end(), cum.begin(), std::plus<double>());
    return cum;
  }
')

これをピュアにRで書くと

cumsum(rnorm(10^5))

とだけ書けばよい。ちなみに、この2つは同じ結果をちゃんと返す。

> set.seed(100)
> randomWalk(10)
 [1] -0.5021924 -0.3706612 -0.4495783  0.4372065  0.5541778  0.8728079  0.2910172  1.0055499  0.1802905
[10] -0.1795716
> set.seed(100)
> cumsum(rnorm(10))
 [1] -0.5021924 -0.3706612 -0.4495783  0.4372065  0.5541778  0.8728079  0.2910172  1.0055499  0.1802905
[10] -0.1795716

計算速度の比較をすると・・・

#結果比較
> library(rbenchmark)
> res <- benchmark(randomWalk(10^5), cumsum(rnorm(10^5)), order="relative")
> res[,1:4]
                 test replications elapsed relative
2 cumsum(rnorm(10^5))          100    1.49    1.000
1    randomWalk(10^5)          100    2.63    1.765

あらら、Rの方が速いよと、1.77倍速いよと。これはRの例の方が、直にCで書かれた部分だけ呼んでるからだろうと推測。
もうちょいいらんメモリ消費を削って、Rcpp版を書きなおして速くすると・・・以下のような感じ。

library(Rcpp)
sourceCpp(code='
  #include<vector>
  #include <functional>
  #include <numeric>
  #include <Rcpp.h>

  using namespace Rcpp;
  // [[Rcpp::export]]
  NumericVector randomWalk(int size)
  {
    RNGScope scope;
    NumericVector rand(Rcpp::rnorm(size));
    std::partial_sum(rand.begin(), rand.end(), rand.begin(), std::plus<double>());
    return rand;
  }
')

で、結果を比べると10%程度Rcpp版の方が速いというレベルか。

> res <- benchmark(randomWalk(10^5), cumsum(rnorm(10^5)), order="relative")
> res[,1:4]
                 test replications elapsed relative
1    randomWalk(10^5)          100    1.32    1.000
2 cumsum(rnorm(10^5))          100    1.49    1.129

std::shared_ptrを噛ませても、親クラスに子クラスのオブジェクト突っ込める

「あれ?これOKだったっけ?」という話の確認で、具体的には「Baseを継承したDerivedのshared_ptrをBaseのshared_ptrに入れれたっけな?」という。
よく考えると、これが出来ないとshared_ptr使って動的多態できなくなるので、出来て当たり前か…

↓サンプルコード、ちゃんと代入出来てる&コンパイル通る&欲しい結果になる。

#include <iostream>
#include <memory>
//適当な継承関係のあるクラス
class Base
{
public:
	void sayHello(){ std::cout << "Hello" << std::endl; }
};
class Derived : public Base{};
//Baseを引数にhelloと言わせる関数
void sayHello(const std::shared_ptr<Base> & b){ b->sayHello(); }
//main
int main()
{
	std::shared_ptr<Base> b(std::make_shared<Base>());
	std::shared_ptr<Derived> d(std::make_shared<Derived>());
	sayHello(b);
	sayHello(d);
	return 0;
}

autoとdecltypeで楽がしたい

C++03の世界では、mapのiterator作ろうとするとこれえらい量を冗長な感じでタイプしなければならなかったが、C++11から導入されたdecltype使うとなんか…C++なのにふつくしい、そんなコードに思えてきたっつー話。
例えば

std::map<int, std::string>::iterator getIterator(std::map<int, std::string> & x)
{
	return x.find(0);
}

って書いていたものが、

auto getIterator11(std::map<int, std::string> & x) -> delete_reference<decltype(x)>::type::iterator
{
	return x.find(0);
}

と書けるわけです。これでmapのstd::stringが突然charに変わるような仕様変更でも関数の引数を直すだけでよいわけです。冗長性の排除ふつくしい
delete_referenceは自作の参照外し構造体。こちらもお勉強兼ねて作っただけで、通常は標準ライブラリのremove_reference使えばよい。

↓コード全体

#include <iostream>
#include <string>
#include <map>

//referenceはずし
template<class T>
struct delete_reference
{
	typedef T type;
};
template<class T>
struct delete_reference<T&>
{
	typedef T type;
};

std::map<int, std::string>::iterator getIterator(std::map<int, std::string> & x)
{
	return x.find(0);
}
auto getIterator11(std::map<int, std::string> & x) -> delete_reference<decltype(x)>::type::iterator
{
	return x.find(0);
}
int main()
{
	std::map<int, std::string> map;
	map[0] = "A";
	map[1] = "B";
	map[2] = "C";
	std::map<int, std::string>::iterator x1 = getIterator(map);
	auto x2 = getIterator11(map);
	std::cout << "x1:" << x1->second << std::endl;
	std::cout << "x2:" << x2->second << std::endl;
	return 0;
}

実行結果は

x1:A
x2:A

どちらも同じだよと。