Rからc++コードを呼び出すRcpp+inlineを試す

従来、RをCで拡張するときは「Rの基礎とプログラミング技法(→amazon)」にもあるとおり、R.h をインクルードして SEXP (S EXPression) 型と格闘するイメージだけど(何度か実践で使ったけど気楽にやろう、ってものではなかった)、Rcpp + inline を使うと魔法のように簡単に R と c++ を連携させることができるようだ。

なんとインラインで c++ を書いて、それを直接実行する、というすごい実装(inlineパッケージを使った場合。Rcpp は R.h のラッパー)。

使い道としてはMCMCとかやたら遅いアルゴリズムを実装するときに便利そう(C++の中でRの豊富な乱数生成が使えたりするから)。

他にもループが遅い時にループの深いところを c++ で書き換えてしまうとか。かなり気楽だからこそできるワザ。

勉強を兼ねて適当なテスト関数を書いてみた。
“Rからc++コードを呼び出すRcpp+inlineを試す”の続きを読む

[c++11] ネストしたクラスでの decltype

c++11 の decltype について。入れ子クラスの場合のちょっとした問題。下のサンプルコードの一番最後の行でコンパイルエラーが出る。

struct Outer {
  int x;
  struct Inner {
    int y;
  };
};


int main(){
  Outer outer;  // OK
  Outer::Inner inner;  // OK

  decltype(outer) outer2;  // Valid for c++11
  decltype(inner) inner2;  // Valid for c++11

  decltype(outer)::Inner inner3; // Error!!!!
}

コンパイル。

$ g++ decltype_test.cpp --std=c++0x
decltype_test.cpp: In function ‘int main()’:
decltype_test.cpp:16: error: expected initializer before ‘inner3’

decltype(outer) は Outer と展開されるので decltype(outer)::Inner は Outer::Inner と展開されて欲しいところだけど、できないみたいですね。g++ のバージョンは 4.4.3 と 4.5.2 でためしてみました。

ほとんどのケースでは auto を併用すれば問題は解決するのだけど、たまに困ります。これは仕様なのか、未実装なのかはわかりませんが。

Eigen ー C++で線形代数を!(3)

前回(といってもずいぶん前のことですが)は eigen の基本的な使い方について書きました。

Eigen ー C++で線形代数を!
Eigen ー C++で線形代数を!(2)

今回は連立方程式を解くための一つの方法である LU 分解を eigen から使う方法について書いてみたいと思います。

毎度のことですが、より詳細な情報についてはマニュアルを読んでください。英語ですがわかりやすく書いてあります。ちなみに今回使っている Eigen のバージョンは 3.0 です。

Eigen を使うと非常に簡単に LU 分解を計算することができます。余計なお世話かもしれませんが、こんなに簡単だとついついブラックボックスとして使いがちです。しかし、LU分解による求解は数値解析の観点からはベストの方法であるとは限りません。特に、疎行列のケースでは(CG法系列の)別の解法を使ったほうがよいでしょう。Eigen も sandbox レベルでは BiCG などの疎行列用の解法が存在しています。Eigen の疎行列ライブラリを使って自分で実装してもそんなに問題はないと思います。

いきなり御託を並べてしまいましたが、気を取り直して、先に進みましょう。
“Eigen ー C++で線形代数を!(3)”の続きを読む

Eigen ー C++で線形代数を!(2)

(このエントリーは旧ブログからの転載です)

ここではC++用線形代数ライブラリ eigen を個人的に習得するために作ったサンプルプログラムを公開しています(一部、本家のチュートリアルと重複しています)。使ってみて感じたのは C++ という言語で ruby や R 的な配列の柔軟さを実現していてすごいなーということでした。ブロードキャストの例に出てくるメソッドチェーンなんて結構感動(式テンプレートのおかげで、おそらくパフォーマンスは落ちない)。すばらしいライブラリです。以下、ライブラリの基本的な説明です。固有値計算などいわゆる線形代数特有の演算についてはいずれ書きます。

【追記】 続きを書きました。Eigen ー C++で線形代数を!(3)

“Eigen ー C++で線形代数を!(2)”の続きを読む

Eigen ー C++で線形代数を!

このエントリーおよび続くエントリーで、C++ で書かれた線形代数ライブラリ Eigen を使うためのメモなどを書いてみようと思います(このエントリーは旧ブログからの転載です)。

とはいっても本家のほうに非常にわかりやすいマニュアルが整備されているので、「この部分はマニュアルを見てね」的な部分も出てくるかもしれませんが。
最終的には微分方程式とか最適化とかの応用問題について書いてみたいと思います。

(現時点(2011/2)で私が触っているのはバージョン 2.0 ですが、もうすぐバージョン 3.0 が出るようです)
“Eigen ー C++で線形代数を!”の続きを読む

【C++】 基本型・クラス間の型変換を定義するには?

c++ では static_cast<double>(int) とかやると型を変換できますが、それを自分で定義したクラスに対して行う方法について書いてみます。
コンストラクタを使ってやる方法もありますが、この方法だと int, double などの基本型(組み込み型)には使えません。たとえばあるクラス X を int 型にキャストしたい場合などですね。

X::operator int() を定義する

こういうときは次のような方法が有効です。

struct X {
  operator int(){ /* 変換のためのコード */}
};

このような operator T を定義しておくと、

X x;
int y=x;
int z=static_cast<int>(x);

というような書き方ができるようになります。

サンプル

上で書いてあることとほとんど同じですが、コンパイルできるコードを一通り書いておきます。
さまざまな型から Hoge クラスにキャストしたり、また、Hoge クラスから別の型へのキャストをするようなサンプルになっています。

#include &lt;iostream&gt;
#define castmsg(from, to) std::cout << "Cast from "#from" to "#to << std::endl;
struct Fuga {
  int fuga;
};
struct Hoge {
  int hoge;
  operator int(){   // Hoge -> int へのキャスト
    castmsg(Hoge,int);
    return hoge;
  }
  operator Fuga(){  // Hoge -> Fuga へのキャスト
    castmsg(Hoge,Fuga);
    Fuga f;
    f.fuga = hoge;
    return f;
  }
  Hoge(int h){      //  int -> Hoge へのキャスト
    castmsg(int,Hoge);
    hoge=h;
  }
  Hoge() : hoge(0) { }  // デフォルトコンストラクタ
};

int main(){
  Hoge hoge;
  hoge.hoge = 1;
  std::cout << "=== Type cast test (by static_cast) ===" << std::endl;
  /* Hoge -> int */
  int h = static_cast<int>(hoge);
  /* Hoge -> Fuga */
  Fuga fuga = static_cast<Fuga>(hoge);
  /* int -> Hoge */
  Hoge hogera = static_cast<Hoge>(4);
  std::cout << "Hoge->int  : " << h << std::endl;
  std::cout << "Hoge->Fuga : " << fuga.fuga << std::endl;
  std::cout << "int->Hoge  : " << hogera.hoge << std::endl;
  std::cout << "=== Type cast test (implicit) ===" << std::endl;
  /* Hoge -> Fuga */
  fuga = hoge;
  /* int -> Hoge */
  hogera = 6;
  std::cout << "hogera: " << hogera << std::endl;    // ここでも暗黙の Hoge->int が使われる
  //std::cout << fuga << std::endl;    // これはエラー!(Fuga->int は存在しない)
  std::cout << "fuga.fuga: " << fuga.fuga << std::endl;
}

/* 出力 */
/*
=== Type cast test (by static_cast) ===
Cast from Hoge to int
Cast from Hoge to Fuga
Cast from int to Hoge
Hoge->int  : 1
Hoge->Fuga : 1
int->Hoge  : 4
=== Type cast test (implicit) ===
Cast from Hoge to Fuga
Cast from int to Hoge
hogera: Cast from Hoge to int
6
fuga.fuga: 1
*/

暗黙の型変換を抑制する

上の例で言うと、ユーザー定義クラスの Hoge に int が = 演算子を使って代入できるようになっていますが、場合によってはこの挙動は危険だ、という場合もあるでしょう。そんなときはコンストラクタ Hoge(int) の前に explict 修飾子をつけましょう。
これをつけることで

Hoge x=1;

のようなコードをコンパイル時エラーにすることができるようになります。

[C++] 列挙体 (enum) から名前(文字列)を取得するマクロ (2)

前回、列挙体 (enum) から名前(文字列)を取得するマクロについて書きました。前回のバージョンは

enum color {RED=1, BLUE=2, GREEN=4};

のような、値の変更に対応していませんでした。というわけで書き換え。

“[C++] 列挙体 (enum) から名前(文字列)を取得するマクロ (2)”の続きを読む

[C++] 列挙体 (enum) から名前(文字列)を取得するマクロ

C/C++ には列挙体があります。列挙体のタグは内部的には整数型として扱われますが、時々、タグの名前を文字列として取得できると便利だな、と思うことがあります。

マクロを使うと割と簡単に実現できます。下のソースコードの DECLARE_ENUM がそのマクロです。可変引数マクロを使っています。その下の main 関数の中では color という列挙体と colorName という文字列配列をDECLARE_ENUM マクロを使って同時に宣言しています。

この手のことはもちろん、文字列配列をまじめに定義すればよいのでどのくらい意味があるかはわかりませんが…
まあ、二度同じことを書かなくてよいのでタイピングは減るということがメリットでしょうか。

#include <iostream>
#include <string>
#include <vector>
#include <boost/foreach.hpp>
#include <boost/algorithm/string.hpp>

using namespace std;

#define DECLARE_ENUM(enumname, enumstr, ...)   \
enum enumname {__VA_ARGS__};                   \
boost::algorithm::split(enumstr, #__VA_ARGS__, boost::is_any_of(",")); \
do { BOOST_FOREACH(std::string &s, enumstr){ boost::trim(s); } } while(0)

int main()
{
  /* color is enum / colorName is vector<string> */
  vector<string> colorName;

  DECLARE_ENUM(color, colorName, RED, BLUE, GREEN);

  cout << RED   << ": " << colorName[RED] << endl;
  cout << BLUE  << ": " << colorName[BLUE] << endl;
  cout << GREEN << ": " << colorName[GREEN] << endl;
}

出力結果

0: RED
1: BLUE
2: GREEN

上の方法の欠点は定義をグローバルな場所に書くことができないことです。列挙体はグローバル変数として定義されることも多いのでこれは問題かもしれません。しかし、すこし書き換えればこの問題は回避できます。

“[C++] 列挙体 (enum) から名前(文字列)を取得するマクロ”の続きを読む

lp_solve で set_presolve したときの目的関数のリセット

以前、lp_solve で set_presolve したときの解の取得というエントリーを書きました。set_presolve は LP を解くための前処理の設定です。これが設定されると solve 関数の実行の際にいくつかの変数が配列から消去されてしまうため、インデックスの対応関係が崩れてしまいます。
前回のエントリーでは solve をした後に、インデックスの対応関係の崩れた lprec オブジェクトからどのように解を取り出せばよいのかについてでした。

今回のエントリーは「lprec オブジェクトのリサイクル」についてです。どういう状況を想定しているかというと、一度 solve が実行された lprec オブジェクトに対して

  • 制約条件はそのままで
  • 目的関数のみを変更する

という、一見すこしマニアックな設定です(実はこの状況は線形制約のシンプルな凸最適化アルゴリズムである、Frank-Wolfe 法を適用する場合に出現します)。

lp_solve ライブラリには set_obj_fn という関数がありますが、前述の理由でこの関数は一度 solve が実行された lprec オブジェクトに対しては使用することができません。

“lp_solve で set_presolve したときの目的関数のリセット”の続きを読む

lp_solve で set_presolve したときの解の取得

lp_solve は C/C++ から使える高性能の線形計画(LP)ソルバーだ。添え字の設定などに多少癖はあるが、中規模程度の問題ならほとんど問題なく利用できるし、なによりフリーなのがうれしい。

この lp_solve には、あらかじめ自明な解や一次従属な制約条件を取り除いてくれるset_presolve という前処理のための関数がある。前処理を行うかどうかは任意(デフォルトではしない)が、することを選択した場合は注意が必要だ。

解は get_variables という関数で取得できるが、前処理されて消去された変数はアルゴリズムの実行時には効率を優先するためにメモリから消されてしまうため、そのままではオリジナルの問題の変数の添え字に正しく対応する値を返してくれない。

“lp_solve で set_presolve したときの解の取得”の続きを読む