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使ってパフォーマンスを測ったらに分類される…?
こういうどうでもいい関数があったとする。
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扱いされる模様。
- Digging into R profiling information - Stack Overflow
- [R] Naming functions for the purpose of profiling
Hadley、パフォーマンス計測用のパッケージ、複数作ってるっぽい。
- GitHub - hadley/lineprof: Visualise line profiling results in R
- GitHub - hadley/profr: An alternative profiling package for R
最終コミットから察するにlineprofを伸ばしていくっぽい。そういえば昔試していた。
呼び出し関係を可視化するproftoolsなんてのもあるのか。
金融商品を記述するためのDSL(Domain Specific Language, ドメイン特化言語)っぽいもんをC++でやってみた(※元ネタ『Composing Contracts: An Adventure in Financial Engineering』)
これは一体何?
金融商品を記述するためのDSL(Domain Specific Language, ドメイン特化言語)を作って、それを実際評価(evaluation)して金融商品のお値段を計算してみましょうという枠組み?のお話。元の本とかだとHaskellなんかで実装するのがセオリーっぽい。自分の呟きでアレですが、Bloombergも既にOCamlで実装した模様。
前々から噂はあったOCamlライクに金融商品を書いて評価できるのブルームバーグのデリバティブライブラリ出てたのか!Excelからも叩ける見たいだしDLIB <GO>するっきゃねぇ! RT BLOOMBERG(PDF) http://t.co/12qXtpG2sK #メモランダム
— 第三使徒テラモナギエル (@teramonagi) 2014, 10月 15
オレオレ実装
割引債とヨーロピアンで力尽きました。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; }
参考
- http://contracts.scheming.org/
- Lexifiの元の論文
- Lexifiの奴のパワポ解説ぽいもの1
- Lexifiの奴のパワポ解説ぽいもの2
- この話を元にしたD論 on GPU実装
- GitHub - channingwalton/scala-contracts: A port of http://web.archive.org/web/20130326233424/http://contracts.scheming.org/- an implementation of Composing contracts: an adventure in financial engineering, by Simon Peyton Jones
- F# Code: Composing Contracts Part 1: Data Types
*1:当記事下部、参考を参照
カジュアルにC++11 in R
以下のプラグイン
plugins=c("cpp11")
をcppFunctionに指定するだけ。Winでもいけた。C++11なんで戻り値をdecltypeにしてやろうと思ったら、できなかった。windowsのgccが古いせいかな?
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
どちらも同じだよと。