Rcpp内で日付(Date)を取り扱う

Rcpp内でDate型を取り扱いたいときにはDateVectorクラスを使用する。
ただ、

によると、Rcpp内で日付を扱うためのDateVectorクラスはそのうちdeprecatedになるらしいが、まぁ、まだあるので使う。

このベクトルの各要素はDateクラスのオブジェクトになっているので、そこは適当にこのドキュメントを見てこいつらを使ってやる。

例1

ためしに「引数として渡したDateのベクトルの日付(8月17日なら17)数分だけ日付を進めた日付ベクトル」をlistとして返すコードを書いてやると、

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List sample1(DateVector x) {
  List result(x.size());
  for(int i = 0; i < x.size(); ++i){
    DateVector date(x[i].getDay());
    for(int j = 0; j < x[i].getDay(); j++){
      date[j] = x[i] + j;
    }
    result[i] = date;
  }
  return result;
}

になる。
実際、動かしてみるとちゃんと動く。

> library("Rcpp")
> sourceCpp("hoge.cpp")
> sample1(x)
[[1]]
[1] "2016-08-28"

[[2]]
[1] "2016-08-29"

[[3]]
[1] "2016-08-30"
> sourceCpp("hoge.cpp")
> x <- Sys.Date() + 1:3
> x
[1] "2016-08-28" "2016-08-29" "2016-08-30"
> sample1(x)
[[1]]
 [1] "2016-08-28" "2016-08-29" "2016-08-30" "2016-08-31" "2016-09-01" "2016-09-02" "2016-09-03" "2016-09-04" "2016-09-05"
[10] "2016-09-06" "2016-09-07" "2016-09-08" "2016-09-09" "2016-09-10" "2016-09-11" "2016-09-12" "2016-09-13" "2016-09-14"
[19] "2016-09-15" "2016-09-16" "2016-09-17" "2016-09-18" "2016-09-19" "2016-09-20" "2016-09-21" "2016-09-22" "2016-09-23"
[28] "2016-09-24"

[[2]]
 [1] "2016-08-29" "2016-08-30" "2016-08-31" "2016-09-01" "2016-09-02" "2016-09-03" "2016-09-04" "2016-09-05" "2016-09-06"
[10] "2016-09-07" "2016-09-08" "2016-09-09" "2016-09-10" "2016-09-11" "2016-09-12" "2016-09-13" "2016-09-14" "2016-09-15"
[19] "2016-09-16" "2016-09-17" "2016-09-18" "2016-09-19" "2016-09-20" "2016-09-21" "2016-09-22" "2016-09-23" "2016-09-24"
[28] "2016-09-25" "2016-09-26"

[[3]]
 [1] "2016-08-30" "2016-08-31" "2016-09-01" "2016-09-02" "2016-09-03" "2016-09-04" "2016-09-05" "2016-09-06" "2016-09-07"
[10] "2016-09-08" "2016-09-09" "2016-09-10" "2016-09-11" "2016-09-12" "2016-09-13" "2016-09-14" "2016-09-15" "2016-09-16"
[19] "2016-09-17" "2016-09-18" "2016-09-19" "2016-09-20" "2016-09-21" "2016-09-22" "2016-09-23" "2016-09-24" "2016-09-25"
[28] "2016-09-26" "2016-09-27" "2016-09-28"

例2

次に、整数値としていったん受け取って、それをDateに直す方法。
これはRcppの中で単にRのas.Dateをimportして呼び出せるようにすればOK。
なんでこんな面倒なことをしないといけないのかというと、Dateクラスのオブジェクトには、数値”ベクトル”に対する加算演算子が定義されていないからだ。

#include <Rcpp.h>
using namespace Rcpp;
Function asDate("as.Date");
// [[Rcpp::export]]
List sample2(IntegerVector x, IntegerVector day) {
  List result(x.size());
  for(int i = 0; i < x.size(); ++i){
    result[i] = asDate(x[i] + seq_len(day[i]), "1970-01-01");
  }
  return result; 
}

これもちゃんと動く。

> sourceCpp("hoge.cpp")
> sample2(x, 1:3)
[[1]]
[1] "2016-08-29"

[[2]]
[1] "2016-08-30" "2016-08-31"

[[3]]
[1] "2016-08-31" "2016-09-01" "2016-09-02"

数値の頭に0を詰めて桁を揃える

C++版をやりたい。

#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>

int main()
{
	// 現状保存
	std::ios::fmtflags curret_flag = std::cout.flags();
	
	//123の頭に5個0を詰めて8桁にする
	std::ostringstream ss;
	ss << std::setw(8) << std::setfill('0') << 123 << "\n";
	std::string s(ss.str());
	std::cout << s;

	// 戻し
	std::cout.flags(curret_flag);
	return 0;
}

実行結果

00000123

C++/C言語の実行時間計測

したい。
関数によってミリ秒だったり秒だったり、OSによって実時間だったりCPU時間だったりする点に注意。

#include <iostream>
#include <vector>
#include <chrono>
#include <time.h>

int main()
{
	const int size = 1000 * 1000;

	//C++11版, ミリ秒単位、実時間ベース
	std::vector<int> x;
	auto chrono_start = std::chrono::system_clock::now();
	for (int i = 0; i < size; ++i) {x.push_back(i);}
	auto chrono_end = std::chrono::system_clock::now();
	std::cout << "Elapsed time(C++):" << std::chrono::duration_cast<std::chrono::milliseconds>(chrono_end - chrono_start).count() << "[ms]" << std::endl;

	//C版(1), CLOCKS_PER_SEC 分の1秒単位、実行時間ベース
	std::vector<int> y;
	clock_t clock_start = clock();
	for (int i = 0; i < size; ++i) { y.push_back(i); }
	std::cout << "Elapsed time(C-1):" << 1000.0*static_cast<double>(clock()- clock_start) / CLOCKS_PER_SEC << "[ms]" << std::endl;

	//C版(2), 秒単位
	time_t time_start, time_end;
	time(&time_start);
	for (int i = 0; i < size; ++i) { y.push_back(i); }
	time(&time_end);
	std::cout << "Elapsed time(C-2):" << 1000*difftime(time_end, time_start) << "[ms]" << std::endl;

	return 0;
}

実行結果
>||cpp|
Elapsed time(C++):1241[ms]
Elapsed time(C-1):1207[ms]
Elapsed time(C-2):1000[ms]
|

数値→文字列をstd::to_string関数でなんとかする

std::to_string関数、引数1個しかないので調整効かないのがあれだったがなんとかする。substrでいらんところを削って、小数点はXに置換しておいた。

#include <iostream>
#include <string>
#include <regex>

int main()
{
	const double pi = 3.1415926;
	std::cout << std::to_string(pi) << std::endl;
	std::cout << std::to_string(pi).substr(0, 1) << std::endl;
	std::cout << std::to_string(pi).substr(0, 2) << std::endl;
	std::cout << std::to_string(pi).substr(0, 3) << std::endl;
	std::cout << std::to_string(pi).substr(0, 4) << std::endl;
	std::cout << std::to_string(pi).substr(0, 9) << std::endl; 
	//拡張子をXに置換
	std::cout << std::regex_replace(std::to_string(pi), std::regex("\\."), "X") << std::endl;
	return 0;
}

結果

3.141593
3
3.
3.1
3.14
3.141593
3X141593

ERROR:invalid initialization of non-const reference of typeを直す

Visual Studio(MSVC)からGCCに移植する際に食らったエラーメモ。

#include <iostream>
template<class T_>
int hoge(int value, T_ & x)
{
	return x(value);
}

struct PlusOne
{
	int operator()(int x){ return x + 1; }
};
int main()
{
	std::cout << "hoge: " << hoge(1, PlusOne()) << std::endl;
	return 0;
}

上の適当コード、これはMSVCコンパイルできるが、GCCだと

$ g++ main.cpp
main.cpp: 関数 ‘int main()’ 内:
main.cpp:15:44: エラー: invalid initialization of non-const reference of type ‘PlusOne&’ from an rvalue of type ‘PlusOne’
  std::cout << "hoge: " << hoge(1, PlusOne()) << std::endl;
                                            ^
main.cpp:4:5: 備考: in passing argument 2 of ‘int hoge(int, T_&) [with T_ = PlusOne]’
 int hoge(int value, T_ & x)
     ^

となりアウト。これを直すにはhoge関数の引数をconst参照にするようにして、かつ、それに合わせてPlusOne構造体のoperator()もconstにすればOK。

#include <iostream>
template<class T_>
int hoge(int value, const T_ & x)
{
	return x(value);
}

struct PlusOne
{
	int operator()(int x) const{ return x + 1; }
};
int main()
{
	std::cout << "hoge: " << hoge(1, PlusOne()) << std::endl;
	return 0;
}

あるいは

PlusOne plus_one;
std::cout << "hoge: " << hoge(1, plus_one) << std::endl;

として一旦実態持たせておくでもOKだった。

今回はちゃんとするなら確かにconst入れるのが筋なので、全部書き直すためにずっとconst, constしていた。

C言語の関数ポインタを取るインターフェイスにRcpp::Functionを渡したい

最近だとtemplateの魔力を使ったC++を使えばこんなことする必要はないですが、古いC言語で書かれたライブラリなんかに、関数ポインタを引数として渡す必要がある場合のお話。

以下のコードでは、関数ポインタ(typedef double(*Simulate)(double x))を引数にとるold_c_function関数にRcpp::Functionを投げ込んでみた例。old_c_function関数自体は適当な処理(x0 + f(x0))を実行する関数としています。

library(Rcpp)
sourceCpp(code='
  #include <Rcpp.h>
  typedef double(*Simulate)(double x);
  //Old-C function
  double old_c_function(double x0, Simulate f)
  {
    return(x0 + f(x0));
  }
  // [[Rcpp::export]]
  double new_interface(double x, Rcpp::Function f)
  {
    return old_c_function(x, &f);
  }
')

これをコンパイルしようとすると、これは無常にもというか、当たり前ではありますが「型変換できねぇよ!」と怒られる。

g++ -m64 -I"C:/R/R-32~1.0/include" -DNDEBUG     -I"C:/R/R-3.2.0/library/Rcpp/include" -I"C:/Users/teramonagi/AppData/Local/Temp/RtmpMR3pkQ"  -I"d:/RCompile/r-compiling/local/local320/include"     -O2 -Wall  -mtune=core2 -c file16514334dad9.cpp -o file16514334dad9.o
file16514334dad9.cpp: In function 'double new_interface(double, Rcpp::Function)':
file16514334dad9.cpp:12:32: error: cannot convert 'Rcpp::Function* {aka Rcpp::Function_Impl<Rcpp::PreserveStorage>*}' to 'Simulate {aka double (*)(double)}' for argument '2' to 'double old_c_function(double, Simulate)'
file16514334dad9.cpp:13:3: warning: control reaches end of non-void function [-Wreturn-type]

これをなんとかしたい、速度とか気にしてRcpp使ってるんだけど、遅くなってもいいから無理やりRcpp::Function使いたいって時は、1枚コンバータ(アダプタ?)を噛ませばいい。staticとしてactionメソッドを取っているのがミソで、普通のメンバー関数だとポインタが変換できなくてダメ。

sourceCpp(code='
  #include <Rcpp.h>
  typedef double(*Simulate)(double x);
  //Old-C function
  double old_c_function(double x0, Simulate f)
  {
    return(x0 + f(x0));
  }
  //Converter for Rcpp::Function <--> Simulate
  struct Converter
  {
    static double action(double x){return Rcpp::as<double>((*f_)(x));}
    static Rcpp::Function* f_;
  };
  Rcpp::Function* Converter::f_ = NULL;
  // [[Rcpp::export]]
  double new_interface(double x, Rcpp::Function f)
  {
    Converter::f_ = &f;
    return old_c_function(x, Converter::action);
  }
')

これはコンパイルできて、ちゃんと動く。以下は「3 + 3^2 + 2 = 14」という計算に対応し、答えもあってる。

> new_interface(3, function(x){x^2+2})
[1] 14

線形代数ライブラリEigenの資料まとめ

C++テンプレートメタプログラミングを活用した線形代数ライブラリであるEigenに関連した資料へのLINKまとめ。

公式ドキュメント

とりあえずここを探す。

特にクイックリファレンスガイド良く見る。

講演資料

特に
Eigen a c++ linear algebra library
が、「Expression Template使ってるぜー」とかの説明もあって良かった。

eigen/Cookbook

短いけどCookbook。うっかりshared_ptrに突っ込んでたので、しょっぱなの「Structures containing Eigen types as members」から勉強になったわ…

Eigen ー C++線形代数を!

flyio氏のまとめ。とりあえず読んでおけ。

行列ライブラリEigenのメモ

暗黒通信団?のメモ(PDF)

調べてた時のツイートまとめ

はじめはArmadillo使おうかなと思っていたが、

ということでEigenにした。
また、より高速なBlazeなんてのもあるらしい。こちらはBoost + Blas必須。


↑上記のLINKのコメ欄でEigen勢から結構突っ込まれてたから、Armadilloの方がいいとは言えないってのにここで気がついた。