Coffee-millなるライブラリを（凄く昔に）書いている。これは主に分子動力学計算の結果を読み書き・解析するためのものなのだが、解析用の部分は別に開発しているのでまあほぼ読み書きのために使っている。今やっていることが落ち着いたらもう少しメンテしたい気持ちはあるが、今回は思い出しがてらそれを開発していた時に何を考えていたかを綴っていこうと思う。
分子動力学の結果のファイルには大抵、粒子の種類と位置座標ベクトルのリストのリストが入っている。分子動力学は古典力学をベースにしているので、ある時刻での粒子の位置（と速度）が決まっており、それが瞬間瞬間と続いていく。なので、ある瞬間の全粒子の位置のリストが一つのスナップショットになり、シミュレーションの時系列データはスナップショットのリストになる。なので、位置のリストのリストになる。特定の粒子だけ追い駆けたい場合も多いので、粒子の種類が一緒に出てくるものもある。
さて、これをどう読めばよいだろうか。ライブラリなら、できるだけ型に意味を持たせたい。一番何も考えずに格納するなら、
double trajectory[3][Nparticle][Nstep];
などになるだろうが、これだと後でどの次元が何を意味していたかわからなくなる。人間の記憶力は一般にさほどよくないし、集中力もない。覚えていてもミスをするし、覚えられなければなおさらだ。
ではどうするか。型だ。型を使う。座標は座標ベクトルとして、「数値の配列」以上の意味を持たせるべきだ。また、スナップショット自体もそうだ。スナップショット型を作るべきだ。そしてそれが集まっているものとしてトラジェクトリができるべきだ。
通常、分子動力学のデータは巨大になるので、スタックに置くべきものではないだろう。ヒープに適切にアロケーションして（ヒープアロケーションは遅いが、確保済み領域にアクセスする速度は殆ど変わらない）用いるべきだが、これを手動でmalloc、realloc、freeまたはnew/deleteするのは良くない。それを管理するクラスが必要だ。そして、それはstd::vectorに提供されている。
ここで注意しておきたいのだが、先ほど言ったとおりヒープアロケーションは遅い。適切なreserveやresizeなしに複数回アロケーションを行ってしまうと、速度が落ちる可能性がある。また、ループ内で一時的に使う短い配列としてstd::vectorを使うと、それをコンパイラが最適化で外に出さなかった場合、非常に遅くなる。しかしこれはどちらかというとその機能が作られた目的を知らないことによって起きたミスで、std::vectorが悪いのではない。ループ内の一時配列でmalloc/freeをする人間はおらず、組み込み配列をスタックに積むだろう（あるいは環境によってはallocaを使ってスタックに積む）。これはそのようなものなのである。std::arrayを使うべきだ。
なので、型に意味を持たせるとすると以下のようになるだろう。
typedef std::array<double, 3> coordinate; typedef std::vector<coordinate> snapshot; typedef std::vector<snapshot> trajectory;
だが、ここで問題が生じる。coordinateとして使える型は非常にたくさんある。もし強力な線形代数ライブラリが使いたいなら、EigenやBoost.uBLAS、blazeなどを使いたいだろう。あるいは、自作の線形代数ライブラリを使いたいかもしれない。もっというと、場合によっては線形代数の手法は必要なく単に座標データがあれば事足りるかもしれず、そのような場合に巨大なライブラリをインストールさせるのは無駄だろう。
なので、線形代数ライブラリがなくとも使えて、また線形代数ライブラリを使う場合は直接そのライブラリが提供する型を返すことができたら理想だ。つまり、coordinate型はtemplateにならざるを得ない。そうすれば、floatで十分だしdoubleを使いたくないという要求にも耐えられる。
template<typename coordT> struct snapshot { using coordinate_type = coordT; std::vector<coordinate_type> positions; }; //または、C++11なら単にこうでもよい template<typename coordT> using snapshot = std::vector<coordT>
何かを読むときは以下のようにする。
const auto traj = read_trajectory<Eigen::Vector3d>("sample.dcd");
そして、デフォルトでベクトル型を用意しておき、関数のデフォルトtemplate引数にしておく(ここでいつの間にかC++11の話になっていたことに気づく。C++98には関数templateのデフォルト引数はない)。すると、普段は単なる関数呼び出しになり、何か特別なことが必要なときにだけtemplate越しに型を指定することができる。
const auto traj1 = read_trajectory("sample.dcd"); //デフォルトベクトル型になる const auto traj2 = read_trajectory<Eigen::Vector3d>("sample.dcd"); //Eigenのベクトル型が直接格納されている
非常に良い。
あとは「ベクトル型」として何が来た時にも対応できるようにするという難しさが残る。普通はoperator[]が定義されているものと想定してよいだろうが、vec.xのようにしてアクセスするライブラリがもしかしたらあるかもしれない。
私の自作ライブラリはそこまで考慮してはいないが、これは例えば一番単純には汎用templateアクセサを用意しておいてユーザーにその特殊化を書かせることで解決できる。低次元ベクトルについてたくさんあるユーザー定義型の差異を吸収するために作られたライブラリとして、Boost.QVMがある。使ってみるのも手だし、その手法を学ぶのもよい。
実際には、単に粒子の位置以上の情報も含まれていることが多いので、単なるcoordinateの配列というわけにはいかないことが多い。
その場合は粒子のattributeをセットにした「あるフォーマットの粒子情報型」というものを作ることになるが、その時もcoordinateを固定するよりは、上のようにする方がよいだろう。
型に情報を持たせる。そして、型をユーザーが変更可能にする。これはtemplateによって両立可能になる。
