C++でターミナルに画像を表示させる

これはC++ Advent Calendar 2013 23日目の記事です。

目的:

sshでリモートに接続して計算を行い、その結果を図として得たい場合、
図をダウンロードするしかないのですが、
いちいちダウンロードして画像ビューアに移ってという過程を経ずに
ターミナル上で画像を表示させたいということが結構あります(私は)。
ターミナルを選びますが、Sixel Graphicsというものを使えばこれは可能です。
これをC++で利用する方法について書きます。

前提:

ターミナルをSixel Graphics対応のものに変えること。
これが一番ハードル高いのですが、
WindowsならRloginが良いと思います。
http://nanno.dip.jp/softlib/man/rlogin/
他の環境だとtanasinnですかね。
http://zuse.jp/tanasinn/usermanual-ja.html
xtermも最近対応してるみたいですが、256色表示させるにはちょっとパッチが必要です。
利用可能かどうかは実際使ってみればわかります。
Netpbmがインストールされていればpngからppmへの変換、ppmからsixel形式への変換ができますので

pngtopnm hoge.png | pnmquant 256 | ppmtosixel

で画像を表示できるはずです。pnmquantは256色以下に減色させるために必要。

大まかな流れ:

1. 画像を読み込む
2. 256色に減色する
3. 画像をSixel形式に変換する
4. 標準出力に吐く
etc) screen/tmux対応

1.は何でもOKです。せっかくC++なのでBoost.GILでも使いましょうか。
2.はBoost.GILではできないようなので、pngquantのライブラリを使わせてもらうことにしました。
3.ここが一番面倒くさいとこですね。実はNetpbmのppmtosixelをパクればいいという話もありますが劣化版を再発明します。
4.吐くだけ。(Boost.GILに合うioを作りたいのが本音なのですが時間足りないので無視します)
etc)実はそのままではscreen/tmux上で使えません。しかし回避策はあります。


1. Boost.GILで画像を読みRGBA形式でのデータへのポインタを得る。

この程度ならGILを使うまでもないのですが、libpngのC言語インターフェイスを使うくらいなら結局GILを使ったほうが簡単です。

#include<iostream>
#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/png_io.hpp>
using namespace boost::gil;

int main() {
    rgba8_image_t src;
    png_read_image("./img.png", src);
    auto raw = interleaved_view_get_raw_data(view(src));
    return 0;
}

2. 256色に減色する

GILで読んだpngデータはrgba8の透過色含む32bitですので、
Sixel形式にするには256色に減色させてやる必要があります。つまりパレットカラーです。

もう256色な時代じゃないので一応説明を加えておくと、256色表示の場合は1pixelの色情報を8bitで表します。
しかし単純に8bitをrgbで表そうとすると、rgbそれぞれが2-3bitしか使えないので、各色2-4階調の色しか表現できないことになります。
そうではなくて256色表示の場合は、使用する色の情報をパレットとして登録しておき、
パレットの1番目が(R,G,B)=(..., ..., ..,)、2番目が(R,G,B)=(..., ..., ...)というように
256色まであらかじめ色情報を作っておきます。
そうすると、1pixelの情報はパレットへのインデックスを保持してればいいので8bitで256色を表現できます。
パレットデータを画像データに含ませる必要があるので、データのヘッダ部分はサイズが増えますが全体としてはサイズが大きく減ります。

256色に減色させる方法は昔から色々あるのですが、一番簡単な方法はR,G,Bをそれぞれ[0,255]の範囲を等分割して、一番近い色を拾っていく手法でしょうか。
あるいはヒストグラムを作って、多く使われる色から順番に選んでいくとか。
どちらも単純な分、それほど再現性は高くないです。
一番よく使用されている手法はmedian cutというアルゴリズムでしょう。
color quantizationとmedian cutでググればすぐ出てくるのでアルゴリズムの詳細は省略します。
時間があればこれも再発明しても良いのですが、1時間程度でサクッとできそうになかったので、何かいいライブラリないかなーと検索してみたらpngquantのライブラリ(libimagequant)がBSDライセンスで公開されていました。
http://pngquant.org/lib/
C言語のライブラリですが、C++のライブラリを探す時間が惜しいのでこれ使います。
一応、C++Compatibleと書かれているbranchがあったのでgit cloneした後cpp branchにしています。
サイトに書かれているように、ダウンロードしてmake -C libすればlibimagequant.aができるので、これを一緒にリンクするだけです。

#include<iostream>
#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/png_io.hpp>
#include"pngquant/lib/libimagequant.h"
using namespace boost::gil;

int main() {
    rgb8_image_t src;
    png_read_image("./img.png", src);
    auto raw = interleaved_view_get_raw_data(view(src));
    // 256色に減色してパレットを取得
    auto attr  = liq_attr_create();
    auto image = liq_image_create_rgba(attr,raw,src.width(),src.height(),0);
    auto result= liq_quantize_image(attr,image);
    auto pal  = liq_get_palette(result);

    // 256色に減色した画像を再構築する
    auto buffer = std::vector<unsigned char>(src.width()*src.height());
    liq_write_remapped_image(result, image, &buffer[0], buffer.size());

    liq_attr_destroy(attr);
    liq_image_destroy(image);
    liq_result_destroy(result);
    return 0;
}

これで例えばi番目の色情報は

auto color = pal->entries[i];
unsigned char r = color.r;

という具合に取得できるようになりました。

3. 画像をSixel形式に変換する

Sixel形式はWikipediaに簡単な説明があります。
http://en.wikipedia.org/wiki/Sixel
さらにこのwikipediaページに貼られている
Chris Chiesa, All About SIXELs, 29 September 1990
に詳しい説明が載っています。これを見ればできそうです。
あと、困ったときは前述のNetpbmのppmtosixelとかのソースを読めばいいんじゃないかと思います。
まあSixel形式自体はそれほど複雑じゃないです。

Sixelという名前ははsix+pixelから来ているようですが、
その名の通り6つのpixelを1セットにしています。
この6つのpixelは縦方向に並びます。
つまり、

0b0000001
0b0000010
0b0000100

という順番に並べられた3つのsixelデータは

x方向->
1 0 0
0 1 0
0 0 1
0 0 0
0 0 0
0 0 0

と表されるわけです。白黒なら1と0で画像を表現できます。
さらにsixelでは0b000001,0b000010,...という数値に63を足し、
それぞれの数値はprintableな範囲に収めています。
sixelデータは6bitで表現されているので、本来、数値としては0-63までです。
この数値に63というオフセットを加える事で、63-126の範囲、
つまりアスキーコード表で?から~までの範囲になり、すべてのデータが文字として視認できるようになります。
例えば、上記の3つのsixelデータは
@AC
という文字列で表現されます。

以上がSixelデータの本質で、あとはいくつかの制御シーケンスとデータ圧縮の話になります。
最低限の項目に絞って解説します。

Sixelモードに入る。

まず最初に指定しなければならない制御コードは
DCS(Device Control Sequence)という文字コードです。
いわゆるエスケープシーケンスというやつで、
7bit表現で"ESP P"
あるいは十進法で144、8進数では220。
文字列として表現するなら、"\0x1bP"です。

さらに、DCSコードに続いてSixelデータの情報とヘッダを出力します。
フォーマットは

DCS p1;p2;p3;q

p1,p2,p3はoptionalです。
・p1はピクセルアスペクト比を示すらしいのですが、互換性のために残されているらしく通常は0としておきます。
・p2は背景色の指定です。どうでもいい感じがしたので説明は省略します。
・p3は水平方向のピクセル間の距離です。

これらの文字コードが標準出力に出されると、Sixelモードに入ります。
この後にカラーパレットの情報を出力しなければならないのですが、
これについては後述します。

Sixelモードを出る。

Sixelモードを抜けるにはST(String Terminator)を標準出力に出します。
STは7bit表現で"ESP \"、8進数で234です。

行頭に戻る。

sixelを出力する位置を行頭に戻すには、'$'コードを出します。
次の行には移りません。
行頭に戻した状態で、sixelデータを出力していくとすでに描かれたデータを上書きしていきます。

次の行に移る。

次のsixel lineの先頭に移ります。
sixelは6pixelを1セットとしますので、6pixelだけ下にラインに移ることになります。

カラーヘッダ

Sixelモードに入った後、カラー情報を出します。HLSでも書けますがRGBでの指定を説明します。

例えば、
#3;2;30;40;50

と書けば、3番目のパレットカラーを(R,G,B)=(30,40,50)とするという意味になります。
注意点として、数値は0-100までのPercentageです。
2番目の数値"2"はRGB表現であることを示しています。

カラーの指定

sixelデータを出している途中で

#5

とか出力すると5番目のパレットカラーに切り替わります。

例:

<ESC>Pq
#0;2;0;0;0#1;2;100;100;0#2;2;0;100;0
#1~~@@vv@@~~@@~~$
#2??}}GG}}??}}??-
#1!14@
<ESC>\

というsixelデータだと

    ..............
    ..**..**..**..
    ..**..**..**..        . = color 1, yellow
    ..******..**..        * = color 2, green
    ..**..**..**..
    ..**..**..**..
    ..............

という画像になります。


4. 標準出力に吐く(コード例)

色々と忙しくまだきちんとしたコード作ってないのですが、とりあえず動くソースコードを出します。

#include<iostream>
#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/png_io.hpp>
#include"pngquant/lib/libimagequant.h"
using namespace boost::gil;

const char DCS='\220';
const char ST ='\234';

template<class Pallete>
int write_header(Pallete* pal) {
     const int maxval = 100;
     const int palmax = 255;
     const int colnum = 255;

     printf("%c",DCS);
     printf( "0;0;8q" );   // 水平方向のグリッドサイズは8
     printf( "\"1;1\n" );  /* set aspect ratio 1:1 */
     for(int i=0;i<colnum;++i) {
          auto color = pal->entries[i];
          const int r = 1.*color.r/palmax*maxval;;
          const int g = 1.*color.g/palmax*maxval;;
          const int b = 1.*color.b/palmax*maxval;;
          printf("#%d;2;%d;%d;%d\n",i,r,g,b);
     }
     return 0;
}

template<class Image>
int write_image(Image& img, const int width, const int height) {
     auto bufsize = width*height;

     const char selector[] = {
          0b000001+63,
          0b000010+63,
          0b000100+63,
          0b001000+63,
          0b010000+63,
          0b100000+63,
     };

     for(int i=0;i<height;++i) {
          for(int j=0;j<width;++j) {
               int index = img[j+i*width];
               printf("#%d%c",index,selector[i%6]);
          }
          printf("$");
          if(i%6==5) {
               printf("-");
          }
          printf("\n");
     }
}

int write_terminator() {
     printf("%c",ST);
     return 0;
}

int main() {
     rgba8_image_t src;
     png_read_and_convert_image("./hoge.png", src);
     auto raw = interleaved_view_get_raw_data(view(src));

     auto attr  = liq_attr_create();
     auto image = liq_image_create_rgba(attr,raw,src.width(),src.height(),0);
     auto result= liq_quantize_image(attr,image);
     auto pal  = liq_get_palette(result);

     auto buffer = std::vector<unsigned char>(src.width()*src.height());
     liq_write_remapped_image(result, image, &buffer[0], buffer.size());

     write_header(pal);
     write_image(buffer,src.width(),src.height());
     write_terminator();

     liq_attr_destroy(attr);
     liq_image_destroy(image);
     liq_result_destroy(result);
     return 0;
}

実行例:

時間がないときのprintfの偉大なる手軽さを知りました。。すごく直したいです。
本当はデータ圧縮(といってもただのRun Length)とかした方がいいのですが、まずはシンプルにしておきます。
screen対応は・・・需要があれば書きます。Twitterでreplyください。tmuxは実は私はしりません。PysixelというPython版の作者が詳しいのでそちらに聞いたほうが早いかも。

CEANによる配列操作

はじめに

これはC++ Advent Calendar 2012の21日目の記事です。
本当はスパコン環境でのC++とかODEintについて書こうかと思っていたのですが、時間がなくなってしまったので予定よりも軽めの内容にします。もしかしたらODEintについて追記するかもしれません。しないかもしれません。

...というか書いてる途中で時間切れになったのでまた後で修正します(申し訳ない)

さてこのCEANというのは何かというと、C/C++ Extensions for Array Notationsのことです。Array Notation。つまり配列の記述方法です。残念ながらC++標準では使えません。Intel Compilerでのみ利用できる拡張機能です。Intel Compiler 13で使えることは確認してます。
拡張機能、しかも文法上の拡張なので使うと他コンパイラとの互換性が無くなります。そのため、この記事で紹介する機能は積極的な使用を推奨するものではありません。

Introduction

C++は実行速度がそこそこあるので数値計算分野にも使えるわけですが、現状として私の環境、特にシミュレーションでは数値計算はほぼFortranです。私も今年C++ほとんど書いてません。というか書けませんでした。スパコンでの利用を考えるとC++はなかなか難しい気もします。速度面での問題もあるのですが、やはりFortran数値計算コードを書きやすい感があります。Fortranといっても昔のじゃないですよ、Fortran95(以降)です。ちょっとこの辺の例も書いてみます。(時間なくてコードチェックしてないのでミスしてたらすみません)
たとえば、3次元での計算コードとかだとほとんどの場合3重ループになるかと思います。

for(int i=0;i<NX;++i) {
  for(int j=0;j<NY;++j) {
    for(int k=0;k<NZ;++k) {
         fx[i][j][k] = (p[i+1][j][k]-p[i-1][j][k])*idx
    }
  }
}

Fortranで書くとこうなります。

forall(i=1:NX-1, j=0:NY, k=0:NZ)
   fx(i,j,k)=(p(i+1,j,k)-p(i-1,j,k))*invdx
end forall

並列化も期待できます。実際はもっと複雑な式になるのでforallを使いましたが、本当はこれだけならさらに簡単に書けます。

fx(1:NX-1, :, :) = (p(2:NX, :, :)-p(0:NX-2, :, :))*invdx

はい。コロンを使った書き方が出てきました。別の例を見てみます。めんどくさいのでC++での例は省略。

a(4:6,4:6) = b(4:6,4:6)  !  配列の一部をコピーします。
S = sum(c(c>0)) !  配列cの内、正の値を持つものの合計
d[0:NX:2]  = 1 ! 偶数番目に1を代入 
e(:) = sin(e(:))  !  e(i) = sin(e(i))をeの全要素に対して適用

こんな感じで、とにかく配列操作が強力です。こういう書き方ができるのでC++よりFortranで書いた方がかなり楽なことが多いです。数値計算のコードはあまり複雑にならないですし(多分)、やはりそれ専用の言語は強いです。例に上げた中でいくつかは、C++でもstd::valarrayとかのライブラリを使えば近いことはできます、sliceとかもちゃんとありますし。でもFortranほどシンプルには書けません。ちなみに、こういうコロンを使ったアクセス方法はPythonのnumpyもできますね。MATLABとかもできた気がします。便利です。C++でも使いたいです。FortranじゃなくてC++を使いたいんです。CEANならできます! しかもCEANはvectorizationもしてくれます。つまりSIMDで並列計算されます。言語拡張だけどちょっと触ってみてもいいんでない?

CEANのつかいかた。

最初に書きましたがIntel Compilerを使います。最近のならデフォルトで有効になってるはずです。追加で必要なコンパイルオプションはありません。
基本的には他言語のNotationと同じです。ただし、指定するのは範囲ではなくて、下限(lower bound)と長さ(length)です。

syntax定義:
 [ :  [: ] ]
             [ :  [: ] ] …
Simple example:
A[:]   // Aの全要素
B[6:2] // B[6]からB[8]まで。指定するのは「長さ」
C[:][5] // 多次元配列も可
D[0:8:2] // 偶数番目
E[:]*F[:] // 対応する要素同士で掛けます。内積とか外積ではない。
G[3:2][3:2] + H[4:2][2][6:2] // 次元やサイズがあってればよい。
I[:][:] + 3  // scalarを足したり掛けたりはOK
Build-in function:

リダクションとかは関数が用意されてます。
__sec_reduce_add(J[:] * K[:]) // 内積
ほかにもいろいろ。

後で書きなおす編

(...ここまで書いて時間なくなってきた。後でまた直そう)

残念ながらはポインタか配列型しかとれないのでstd::vectorとかは無理っぽい

//  std::vector<int> A{1,2,3,4,5,6,7,8,9}; // 無理だった
//  std::vector<int> A{9}; // 通るけど A(9)になった。えー
    std::vector<int> A(9); // サイズを指定するときは丸カッコですよね
    int idx[9];
    idx[:] = __sec_implicit_index(0); // idx[i] = i になる。便利。
                                      // 0から始まるという意味じゃなくて0番目のランクという意味です
                                      // 多次元配列だとさらに便利
//  A[1:4] += 1; // NG
//  A[idx] = idx[:]; // NG
    A[idx[:]] = idx[:]; //OK
    for(auto v: A) std::cout << v << ","; //0,1,2,3,4,5,6,7,8

Fortranの例で出した

e(:) = sin(e(:))

はBuiltin関数としてsinが使えるので

e[:] = sin(e[:])

と書けます。Builtinではない任意の関数で同様のことをするには

__declspec(vector) double sincos(double a, double b,){
   return sin(a)+cos(b);
}
//...
A[:] = sincos(B[:], C[:]);

でできます。都合により私の時間がないので、後で書きなおすことにして次のif文で締めます。

CEAN with if statement

いわゆるmaskingができます。流体解析の風上差分とかでは速度の符号によって計算式が変わったりすることもあるので、そういうときに便利かな?

    // CEAN with if statement
    int mask[8] = { 0,1,1,1,0,1,1,0 };
    if(mask[:]==0) { mask[:]=4; }                                                          
    else { mask[:]=8; }
    for(auto v: mask) cout << v << ","; // 4,8,8,8,4,8,8,4,

今回説明したようなことは次のリンク(PDF)にほとんど書いてありますので、
ぶっちゃけこのPDF見たほうが早いですね。時間ができたらまた続き書きます。ODEintは...年内は厳しいかなあ...
Intel Documentation
http://software.intel.com/sites/default/files/m/0/e/c/3/d/28991-Effective_Data_Parallelism_Using_CEAN.pdf

次は@hgodaiさんの「ゲームプログラミングの都市伝説」です。

PythonからFortranのサブルーチンを呼ぶ。

本題の前に

どうやら日記放置しすぎたようですが、それはまあ気にしない。
今回のお題はPythonです。
最近Science with PythonというコミュニティがGoogle groupにできましてですね、Twitterでのハッシュタグは#sciwpyなんですが、せっかく参加したので科学技術計算で残念ながら現役バリバリのFortranPythonから利用する方法を書いてみよーかと。

PythonからFortranを使う利点

Pythonにもnumpyとかscipyとか便利な数値計算のためのライブラリがあるので、軽い数値計算ならPythonだけでも十分いけます。PyPyとかCythonという選択肢もありますし。
ところが、かなり重い計算とか、すでにFortranでコードがある状態だとPythonで書くor書き直す気にはならないのですよ。なぜか?

  • Pythonで速く書けても、CやFortranならもっと速くなるのではないかという疑念が憑き纏う
  • Fortranで書かれたコードがあるならそれ使えばいい。コードがどれだけ良くなっても結果が変わらないなら成果にならない。

だからヘビーな計算だと、Pythonのみってのはやる気にならない。やるならアルゴリズムのテストとかちょっとした実験とか。
とは言ってもですよ、Fortranとか色々不便なのでPythonを使いたいわけですよ。手軽だし。というわけで、

  • 既存のFortranコードをライブラリ化してPythonから利用する
  • 処理の重い部分だけをFortranで書く

という方針がたちます、これによって、

  • 今までのコードを大きく書きなおさずに利用可能
  • Fortranコードのプリプロセス/ポストプロセス(例えば、計算のバッチ処理化とか計算結果のグラフ出力とかいろいろ)をPythonで手軽に記述でウマー

というメリットがあるわけですね!

どうやって呼び出すのか?

過去にはf2pyとかいうものもあったみたいですが、なんだか更新されてないようです。これについては今回はスルーして、Pythonの標準ライブラリctypesを使います。あと、以下コンパイラ依存が強くなるので話をgcc/gfortranに限定します。今回はPython2.7でやります。
ctypesはPython2.5から追加された、C言語の関数を呼び出すためのライブラリです。C言語Fortranを組み合わせたプログラムを書いたことのある人は知ってると思いますが、実はC言語からFortranの関数(サブルーチン)を呼び出すのは結構簡単です。C言語からFortranのサブルーチンを呼ぶときの注意点は

  • サブルーチン名の末尾にアンスコ"_"を付けたものがC言語から呼び出す時の関数名となる。
  • サブルーチンの返値はvoid型
  • 引数は参照(ポインタ)渡し

といったところです。C言語Fortranは簡単に連携でき、C言語Pythonは標準ライブラリctypesで連携できる!
ならば、C言語Fortranのラッパーを書いてやればいける!



と思うかもしれませんが(最初私はそう思いましたが)、C言語でラッパーを書く必要はありません。ctypesでFortran直接呼べます。少なくともgcc/gfortranであれば。

Example 1 基本編

具体的に、次のようなFortranの関数をPythonから呼び出すことを考えます。

! add.f90
subroutine add(a,b)
    implicit none
    integer(8),intent(in) :: a
    integer(8),intent(inout):: b
    c = a + b
end subroutine

やってることは簡単です。bにaを加えるだけです。intentとか無くても動きますがつけます。注意事項ですが、integerのサイズ(32bitか64bitか)が環境依存なので、きちんと指定したほうが良いかと思います。この記事では整数・浮動小数いずれも64bitです。このソースを共有ライブラリとしてビルドします。

$ gfortran -shared -fPIC -o add.so add.f90

これはLinuxの場合の例です。Windowsの場合はfPICを消して、soではなくdllにして下さい。

次にPythonでのコードです。

#!/usr/bin/env python 
#vim:fileencoding=utf8  
from ctypes import * 

addmodule = cdll.LoadLibrary("add.so")

# 引数の型を指定する。今回はintのポインタ型
# Fortranのサブルーチン名"add"にアンスコ"_"をつける(gccでは)
addmodule.add_.argtypes = [ POINTER(c_int64),POINTER(c_int64)]
      
# 戻り値の型を指定する。Fortranのサブルーチンはvoidしか返せない
addmodule.add_.restype = c_void_p  

# 呼び出しに使う引数はctypesの型でラップする
a,b = 10,8 
a = c_int64(a)
b = c_int64(b)
 
# byrefでポインタにして渡す         
addmodule.add_(byref(a),byref(b)) 
print b.value # 18 

詳しい説明はソースコード中に書いたので省きます。実はC言語(というかポインタ)が不慣れなウチの後輩がハマっていたので、それを紹介するためにa,bの初期化で若干回りくどい書き方をしています。というのは、この例ではc_intでラップした変数をbyrefを利用して参照で渡していますが、もし

a,b = 10,8  
addmodule.add_(byref(c_int64(a)),byref(c_int64(b)))
print b  

と呼び出した場合はどうなるか? これ結果は8になります。わかっている人には何でもないことですが、c_int64でラップすると別のオブジェクトが作られるので、その参照をとって渡してしまうと意図した結果が得られません。注意。
あと、よくわかりませんが実はbyref使わなくても動いてしまいました。謎。
byrefではなくpointerを使う方法もありますが、速さで言えばbyrefの方が速いそうです。詳しくはドキュメントへGoってことで。

とまあ、こんな感じで直接Fortanのサブルーチンを呼び出すことができます。
まとめると、

  1. Fortranサブルーチンを共有ライブラリとしてコンパイル
  2. Pythonではライブラリのロード
  3. 関数プロトタイプ(引数、返り値の型)を指定
  4. ctypesのデータ型でラップして呼び出し。Fortranは参照渡しな点に注意。

です!

Example 2 実践編

Fortranのサブルーチンを呼べるようになったわけですが、実際のところ、int型の変数1つとか渡さないわけですよ、数値計算やってる人は。渡すのは主に配列。しかも、Pythonとの連携をするならPython側はnumpyの配列使ってることでしょう。だから、numpyの配列を渡すとこまでやらないといけませんよね?

というわけで、基本編のnumpy配列版。まずはFortranソースコードから

! add_np.f90
#define MAXSIZE 1024 
subroutine add_array(a,b,N)
    implicit none
    integer(8),intent(in) :: N  
    real(8),dimension(0:MAXSIZE),intent(in) :: a
    real(8),dimension(0:MAXSIZE),intent(inout):: b
#ifndef NDEBUG
    if (N>MAXSIZE) then
        print *,"MAXSIZE is too small."
    endif 
#endif 
    b(0:N-1) = b(0:N-1) + a(0:N-1)
end subroutine            

Fortranではほとんど使われない(というかFortranじゃない)プリプロセッサを使って配列のサイズチェックを行なってます。gfortranではコンパイルオプションに-cppを追加することで、C言語プリプロセッサが使えるようになります。
でもBOOST.PPは使えなかったはず...

なんでこんなことするのかというと、Pythonから渡された配列は型としてサイズがわからないからですね。こんな風に大きめに配列をとって、もし小さければ警告を出すなどの工夫が必要になります。プリプロセッサならデバッグが終わったら外せるので便利です。ただし、コンパイラ依存?
渡された配列サイズがわからないので、引数でサイズをとるようにしています。C言語と同じです。
この辺り、もっとマシな方法があればいいのですが・・誰か教えて下さい。

次にPythonの方ですね。numpyにはctypeslibがあるのでそれを使うことにします。使い方はほとんどctypesと同じです。

#!/usr/bin/env python
#vim:fileencoding=utf8
from ctypes import *
import numpy as np

add_np = np.ctypeslib.load_library("add_np.so",".")
add_np.add_array_.argtypes = [
               np.ctypeslib.ndpointer(dtype=np.float64),
               np.ctypeslib.ndpointer(dtype=np.float64),
               POINTER(c_int64),]
add_np.add_array_.restype = c_void_p   
                    
a = np.arange(0.,10.,dtype=np.float64) # 0,1,2,3,4,5,..,9 
b = a*2                                # 0,2,4,6,8,...,18
size = byref(c_int64(b.size)) 
add_np.add_array_(a,b,size) # ndarrayはそのまま渡して良い 
print b # [ 0 3 6 9 12 .. 27 ] 

ほぼExample1と同じですね。numpyの配列は大抵ndarrayとして扱うことになりますので、Fortranにはそのポインタを渡してやります。それがnp.ctypeslib.ndpointer(dtype=np.float64)の部分。私はかなり変数の型に注意を払ってます。
コメントにも書きましたが、ndarrayな配列はそのままFortranのサブルーチンに渡してもOK。
整数型はNG。面倒な場合は整数もndarrayで渡してしまうという手も・・。


numpy配列の場合もあまり説明することないですね。(久しぶりに書いたから疲れたわけではないですよ決して)
こんな感じで、Python+numpy/scipy+Fortran/C,C++が簡単にできるのでScienceな分野だと結構Python強いんじゃないかなーという気がします。

みんなでScience with Pythonしましょう!

あけまして

あけましておめでとうございます。
今年は辰年ということなので、dragonですよね。
というわけでdragonegg(http://dragonegg.llvm.org/)を使えるようにしてみました。

dragoneggとは

dragoneggはGCCプラグイン機能を使ってLLVMによる最適化をしようというものらしいです。
現状、Linuxdarwin(Mac)でしか使えなさそう。

今のところそれほど良いことは無さそうな気がするので
私は使えるようにしただけで終わりました。
インストール方法だけメモっときます。Macのみです。

インストール

gccllvmとdragoneggが必要です。全部macportsで入れます。
落とし穴ですが、macportsgccはltoが有効になってないのでportfileをいじってltoを有効にできるようにします。
本当は野良portfileを作ってゴニョゴニョするのがいいのしれませんが、面倒くさいのでportfileを直接いじります。

$ sudo port edit gcc46

でgcc4.6のportfileを開きます。そんで、

variant lto description {enable the Link Time Optimization} {
    configure.args-append   --enable-lto
}

という記述を追加します。これでmacportsでltoを指定できるようになります。
そしてgcc46をインストール

$ sudo port install gcc46 +lto

うまくインストールできることを願って下さい。次にllvmとdragonegg

$ sudo port install llvm-3.0
$ sudo port install dragonegg-3.0

同時にやってしまってもいいかもしれません。

うまくいけば、dragonegg-LLVMバージョン-xxxとかいうものが使えるようになるはず。
g++なら

$ dragonegg-3.0-g++ -std=c++0x main.cpp -o main

というかんじで。実際はg++に-fplugin=...とかいうオプションつけてるだけなはずなので、g++のコンパイルオプションとかも使えます。


どの程度最適化されるんでしょうね・・・



そうそう、年賀状代わりにdragoneggのロゴ置いときます。

ヌーの人とは関係ないですよ・・

多重forループを書きやすくしたい

のでちょっと書いてみた。C++11限定でBoostが必要。
こんな風に使う。

#include <iostream>
#include "feux/multiple_for.hpp"

typedef std::pair<int,int> Range;

int main() {
    const Range range1 = std::make_pair(0,3);
    const Range range2 = std::make_pair(0,3);
    const Range range3 = std::make_pair(0,3);
    const Range range4 = std::make_pair(0,3);

    // 4重ループ
    feux::multiple_for(
            [](int i,int j,int k,int l){ printf("%d %d %d %d\n",i,j,k,l); },
            range1, range2, range3, range4);
    return 0;
}

16重まではイケるはず
まあお遊びみたいなものなので、使い勝手がいいかは知らない。
ヘッダはこちら

/*
 * 多重forループ文
 */
#pragma once
#include <algorithm>
#include <boost/mpl/int.hpp>
#include <boost/preprocessor/repetition/enum_params.hpp>
#include <boost/preprocessor/repetition/repeat_from_to.hpp>
#include <boost/preprocessor/arithmetic/dec.hpp>
    
namespace feux {
            
template<class Function,class HR, class ... Ranges>
struct multiple_for_holder {
    static void invoke(boost::mpl::int_<1>,const Function func,const HR hr) {
        for(int i=hr.first;i<hr.second;++i) func(i);
    }

#define GEN_FUNCTION(z,n,data) \
    static void invoke(boost::mpl::int_<n>,const Function func,const HR hr,const Ranges ... tr) { \
        auto curried = [func](const int i) { \
            return [i,func](BOOST_PP_ENUM_PARAMS(BOOST_PP_DEC(n),const int I)){ \
                func(i,BOOST_PP_ENUM_PARAMS(BOOST_PP_DEC(n),I)); \
            }; \
        }; \
        for(int i=hr.first;i<hr.second;++i) multiple_for(curried(i),tr...); \
    }
    BOOST_PP_REPEAT_FROM_TO(2,16,GEN_FUNCTION,data)
#undef GEN_FUNCTION
};

template<class Function,class HR, class ... TR>
void multiple_for(const Function func,const HR range, const TR ... ranges) {
    multiple_for_holder<Function,HR,TR...>::invoke(
        boost::mpl::int_<sizeof...(ranges)+1>(),func,range, ranges...);
}
};

本当はfunctorとrangeの引数の順番を逆にしたいのですがそれはまた今度で。
いろいろ試行錯誤しながらだったので、ムダもありそう。

今日の「それC++11/Boostでできるよ!」

http://blog.livedoor.jp/dankogai/archives/51761113.html経由で
http://hirobanex.net/article/2011/12/1324792864を知ったわけですが、
こういうのって需要多いんですね。


これC++11なら標準ライブラリにあります。
C++11にまだ移行できない場合はBoostにあるのでそちらで。

こんな感じで使えます。

#include<iostream>
#include<vector>
#include<random>
int main() {
    std::mt19937 rng;
    std::discrete_distribution<> dist{ 20, 10, 10, 30, 30};
    std::vector<std::string> items{"foo","bar","baz","hoge","moge"};
    std::cout << items[dist(rng)] << std::endl;
    return 0;
}

C++11/boostのrandomは乱数生成器と分布を組み合わせて使えるので、
ものすごく便利です。分布はそのままで、乱数生成アルゴリズムだけ切り替えるとかも簡単にできますし、逆に分布だけ別のものに切り替えることも可。
今回の例のような、重み付き乱択っていうんですか?こういう乱数は結局分布が違うだけなので、分布を作ってやれば十分ということ。
あ、これC++11 Advent Calendarのネタにすればよかったorz


ついでに、本当にこれで出来てるのか確かめてみる。
検証コード:

#include<iostream>
#include<vector>
#include<random>

int main() {
    std::mt19937 rng;
    std::discrete_distribution<> dist{ 20, 10, 10, 30, 30};
    std::vector<std::string> items{"foo ","bar ","baz ","hoge","moge"};
    const int N = dist.probabilities().size();
    const int total=1e6;
    std::vector<int> counts(N);

    for(int i=0;i<total;++i)
        ++counts[dist(rng)];

    for(int i=0;i<N;++i)
        std::cout << "item: " << items[i]
                  << ",probability: " << dist.probabilities().at(i)
                  << ", generated: " << 1.*counts[i]/total
                  << std::endl;
}

結果:

item: foo ,probability: 0.2, generated: 0.199779
item: bar ,probability: 0.1, generated: 0.100002
item: baz ,probability: 0.1, generated: 0.100087
item: hoge,probability: 0.3, generated: 0.300024
item: moge,probability: 0.3, generated: 0.300108

1e6回も回してこの精度がいいか悪いかは知りませんが
ちゃんと出来てますね。