np.arange() in c++ - 忘食忘憂

久しぶりにPythonからc++に戻ったら、 np.arange()やnp.linspace()、そしてnp.append()が使えないことに強烈なもどかしさを覚えた。つまり、vector<double>（あるいはvalarrayかただの配列）の等差数列を始値と差分あたりを指定して生成したり、それらを一発で結合する関数がほしい。

例として、関数double func( double )の値を等間隔で表示したいとき、もしc++にPythonのnp.arange()に対応する関数ARange()があれば以下のように書ける。

  for( auto x : ARange( 10. , 500. , 10. ) )
    cout << x << " "
         << func( x ) << " "
         << endl;

これだけなら普通のfor文でもあまり変わらないが、関数を非等間隔なグリッド上で出力しようとすると通常のfor文だけでは対応が難しくなってくる。そのような需要が、実は数値計算ではかなりある。特に関数の処理に時間がかかるときは、関数の評価を上のように等間隔で行うのは時間の無駄で、関数形に合わせて非等間隔に調整したグリッド上でのみ関数の評価を行うことで計算コストが大幅に削減できる。 Pythonではこのようなグリッドをnp.arange()とnp.append()の組み合わせで簡易的に作っていたので、 c++で同じことをしようとして、はたと手が止まってしまったのである。

それで、c++におけるnp.arange()やnp.append()の代用を探してみたのだが、標準c++ライブラリの中にはそのような役割を果たす関数はどうもなさそうである。

さらに調べていくと、c++上でnumpyライクな書式を実現するライブラリ"NumCpp"
https://dpilger26.github.io/NumCpp/doxygen/html/index.html
なるものが見つかったりして興味津々なのだが、今回は等差数列の生成と結合を行いたいだけなので、このライブラリは大げさすぎるように思う。

こういうとき、どうするのがスマートなのだろうか。

いろいろ考えた末に、自分で実装することにして、ひとまず以下のようなものを書いた。

vector<double> ARange( double bgn , double end , double dif = 1. )
{
  vector<double> tmp(1,bgn);
  for( double val=bgn+dif ; val<end ; val+=dif )
    tmp.push_back(val);
  return tmp;
}

vector<double> LinSpace( double bgn , double end , int div )
{
  vector<double> tmp(div+1);
  double val = bgn;
  double dif = ( end-bgn ) / (double)div;
  for( auto& t : tmp ){
    t = val;
    val += dif;
  }
  return tmp;
}

vector<double> Append( const vector<double>& v1 ,
                       const vector<double>& v2 ,
                       const vector<double>& v3 = vector<double>() ,
                       const vector<double>& v4 = vector<double>() ,
                       const vector<double>& v5 = vector<double>() )
{
  vector<double> tmp = v1;
  tmp.insert( tmp.end() , v2.begin() , v2.end() );
  tmp.insert( tmp.end() , v3.begin() , v3.end() );
  tmp.insert( tmp.end() , v4.begin() , v4.end() );
  tmp.insert( tmp.end() , v5.begin() , v5.end() );
  return tmp;
}

この実装だと、ARangeで差分difを負にできないとか、 Append()の仕様がnp.append()と違うし同時に5つまでしか結合できないとか、もっと効率の良い実装があるかもしれないとかの問題があるのだが、これで冒頭のコードは実現できるし、差し当たっての目的は達成できるので目をつぶろう。

早速、執筆中の論文の数値計算用に書いた関数RateM()の値をファイル出力するのに使ってみる。

  ofstream file( filename );
  auto grid = Append( ARange( 1. , 10. , 1.5 ) ,
                      ARange( 10. , 40. , 5. ) ,
                      ARange( 40. , 100. , 10. ) ,
                      ARange( 100. , 200. , 20. ) , 
                      ARange( 200. , 500. , 40. ) );

  for( auto M : grid )
    file << M << " "
         << RateM( M , C ) << " "
         << endl;

うん、やりたかったことができて大変満足。

さて、こんな感じの関数たちはヘッダファイルにまとめて汎用的に利用するに値するだろうか。またその場合、関数の命名規約やファイル名はどうするべきだろうか。さらに、やるならやるで、もっと柔軟な操作も含めて徹底的にやりたくなるのだが、それなら上記のNumCppを使えという話になるので、そのあたりのさじ加減はどうするのだろうか。

楽しい妄想は続く。