Fortranのバイナリ出力データ(unformatted)をPythonで読む

Fortranで出力したunformatted形式のデータをPythonで読む方法です。 scipyにFortranFileというものがあるので、それを使ってもできます。

scipy.io.FortranFile — SciPy v0.15.0.dev Reference Guide

でもこれソース読めばわかりますが、numpyのfromfileって関数を利用してるだけです。

numpy.fromfile — NumPy v1.8 Manual

fromfileで直接扱ったほうが応用も効くのでこれでやります。環境依存な部分もあるかもしれませんので注意してください。

今回読み込むデータは次のようなFortranソースコードで出力されたとします。

program main
    integer :: a, i,j
    real(8) :: b, arr(4,4)
    a = 4
    b = 3.14
    forall(i=1:4,j=1:4)
        arr(i,j) = 10*i+j
    end forall
    open(20, file="output.dat", form="unformatted")
    write(20) a, b, arr
    close(20)
end program

aは32bit整数、bは64bit浮動小数点数、arrは4x4の2次元配列です。こうすると、output.datというファイルには

という感じでデータが並びます。データの先頭と末尾にwrite文で出力したデータサイズ、今回の場合は

変数aのサイズ(4)+bのサイズ(8)+arrのサイズ(4x4x8) = 140バイト

がSize information(正式な呼び方は知りません)と書いたところに入ります。この領域自体は4バイトです。ファイルに追記していく場合でもwrite分のたびに、サイズが出力されます。

Python側ですが、fromfileのドキュメントを読めば予想がつくかもしれません。こんなかんじになります。

import numpy as np

head = ("head","<i")
tail = ("tail","<i")

dt = np.dtype([head, ("a","<i"), ("b","<d"), ("arr","<16d"), tail])
fd = open("output.dat","r")
chunk = np.fromfile(fd, dtype=dt, count=1)
print "a:",chunk[0]["a"]
print "b:",chunk[0]["b"]
arr = chunk[0]["arr"].reshape((4,4),order="F")
print "arr:",arr

dtypeでデータ構造を指定しています。"<"はリトルエンディアンです。ビッグエンディアン環境で出力した場合は代わりに">"を使います。つまり

dtype	意味
<i	リトルエンディアンの整数
<d	リトルエンディアンの浮動小数点数
<16d	リトルエンディアンの浮動小数点数16個

という意味です。2次元配列とかの指定方法はないので、配列の総数を指定します。データを読み込んだ後、reshapeを使って2次元配列に整形します。order="F"がポイントです。いわゆるC言語系とFortranでは多次元配列の並び方が違うので。

fromfileでcount=1としているのは読み込みをループしないという意味です。今回はwrite文が1回しかないので書かなくても同じですが。複数回write文を呼び出した場合はcountを読み込みたい回数に変更します。chunk[1]に2回目のwrite文で出力したデータが入ります。

これでFortranのデータは読めるようになったので、Pythonでデータを解析してmatplotlibで可視化するとか簡単にできますね！

Haxeを使ってみる

インストール

私はArch Linuxなので

$ yaourt -S haxe

終わり。

C++ Console Application

私はC++erなので、まずC++から。OfficialのGetting started with Haxe/C++通りに進める。 まずC++でやるには、hxcppというものをインストールしなければならないらしい。

$ haxelib install hxcpp

ところが、

> This is the first time you are runing haxelib. Please run haxelib setup first

初回はセットアップしろと。そうですか。

$ haxelib setup

書き込み権限のあるパスを指定しろと聞かれる。デフォルトの場所、ユーザ権限じゃ書き込めないのでとりあえずカレントディレクトリでも指定しとく。場所はあとで考えよう。セットアップが終わったら今度こそ

$ haxelib install hxcpp

hxcppというディレクトリになんか色々入った。48Mくらいある。

とりあえずHello, Worldしましょう。Test.hxというファイルを作って。

class Test {
    static function main() {
        trace("Hello World!");
    }
}

と書く。コンパイル用の設定ファイルも必要らしい。compile.hxmlというファイルを作って

-cpp cpp
-debug
-main Test

と書く。んで、いざコンパイル。

$ haxe compile.hxml

なんかStandard Libraryが見つからんというエラー出たので、ターミナルを再起動してリトライ。とりあえずコンパイラは動いた模様。しかし、

/usr/include/gnu/stubs.h:7:27: 致命的エラー: gnu/stubs-32.h: そのようなファイルやディレクトリはあり ません

なんですかそれは。32bit用のファイルが足りてないようだ。32bitのファイルをインストールしてもいいけど、64bitで行きましょう。compile.hxmlを次のように修正。

-cpp cpp
-debug
-main Test
-D HXCPP_M64

再度

$ haxe compile.hxml

通った。よしよし。cppというディレクトリが作られて実行バイナリがある。動かしてみる。

$ ./cpp/Test-debug

>Test.hx:3: Hello World!

Hello, Worldはサクサク進む。

MatplotlibでScientific notationを使った時のcolorbarの上の部分の文字をカスタマイズ

Matplotlibの目盛操作はFormatterを使ってカスタマイズできます。中でもScalarFormatterというものにはscientific notationという機能があります。これはa x 10^bという表記だと目盛にはaを表示させて10^bは端の方に書いておくというものです。言葉で説明すると面倒くさいので図を出します。

デフォルトで勝手にそうなってるかもしれませんが、これはScalarFormatterのset_scientificでON/OFFを切り替えられます。 例えば上の図を出したソースはこれです。、

import numpy as np
import pylab as pl

x = np.linspace(0,1)
y = np.linspace(0,2)
x,y = np.meshgrid(x,y)
z = 1e-8*(x**2 + y**2)

pl.contourf(x,y,z)
cbar = pl.colorbar()
cbar.ax.tick_params(labelsize=12)
cbar.formatter.set_scientific(True)
pl.show()

上で出した図中の青い枠線で囲った1e-8という文字を大きくしたり、カスタマイズしたいのです。次のように修正することで変更可能です。

import numpy as np
import pylab as pl

x = np.linspace(0,1)
y = np.linspace(0,2)
x,y = np.meshgrid(x,y)
z = 1e-8*(x**2 + y**2)

pl.contourf(x,y,z)
cbar = pl.colorbar()
cbar.ax.tick_params(labelsize=12)
cbar.formatter.set_scientific(True)
text = cbar.ax.yaxis.get_offset_text() # offset_textが操作対象
text.set_fontsize(28)
text.set_color("blue")
pl.show()

こういうの名前がわからないので調べにくい。

Boost.勉強会 #15 札幌を開催しました。

札幌C++勉強会/Sapporo.cpp主催で札幌で２度目となるBoost.勉強会を5/24(土)に開催しました。 運営の裏側の苦労も記録に残しときます。

そもそもBoost.勉強会を開催しようという話はSapporo.cppのSkype運営会議でhotwartermorningさん(id:heisseswasser)の提案で出てきました。最初に話が上がったのはいつだったかな。寒い時期がいいか、暖かくなってきた時期がいいか議題にあがっていて寒いとid:faith_and_braveさんが来ないんじゃね？みたいな話をしてたので結構前だった気がしますね。正直私は色々忙しいことはわかってたのでかなり不安だったのですが、やる気があるならやりましょうかということで承諾しました。ええ当然予想通り今かなり反動で時間がなくて大変ですよ。。ブログ書いてる時間が惜しい。

Sapporo.cppとしても北大を会場に使ったのは初めてです。他にも選択肢はあって２年前のBoost.勉強会 札幌の会場を利用する案もありました。北大を選んだのはネットワーク、プロジェクタが使え、交通アクセスが良く、入場料をとらなければ無料で利用でき、そこそこ広い部屋を借りられたからです。50人以上入れてこのような条件を満たす会場は札幌にはほとんど無いと思います。利用には北大の事務に提出する企画書を作ったりしなければならなかったので、はっきりいって面倒でした。主催しているSapporo.cppは継続的に勉強会を続けていましたし、Boost.勉強会も開催歴がそこそこあり、しかも他大学での開催例があったので企画書が書きやすいのが救いでしたね。江添さんに私がTwitterでお誘いしたことが氏のブログとかスライドに載ってますが、"そこそこ知名度のある人が来るのでぜひ北大で開催したい"と企画書に載せるための伏線だったりします。企画書に載せるには正式な募集より先に発表者を何人か確定させる必要があったので。この辺の作戦は他の運営メンバーも知らないはずですね。

無事企画が通って北大を使えることになったのですが、無料で利用するには学生主体で運営しなければならないという条件を課せられています。私が単独で主催するのがベストだったのですが、Sapporo.cppの顔はhotwartermorningさんですし、私が自由に動けるとも全く思ってなかったので今回だけは主催者を私とhotwartermorningさんの連名という形にさせてもらいました。ただし運営は学生主体です。ここ重要です。企画書作ったりしたのは私ですし実際それなりに働いてますので問題ない。

北大を会場にしてどうだったかと言いますと、メリットもあり、デメリットもありといったところで判断は難しいです。メリットは上に書いたようなこと。デメリットは物品販売ができないとか、ネットワークの問題が発生しやすいことですね。今回、事前にインターネット利用のゲストアカウント申請を行っていて、希望者は北大のネットワークを利用できるはずでした。大変申し訳無かったのですが、申請手続にミスがあったらしく利用できない結果となってしまいました。参加者の皆様には本当に申し訳ないです。不幸中の幸いは会場のクラーク会館は公道から近かったことですね。北大の敷地は広いので、もし工学部とかだったら公道から遠すぎてモバイルルータとかも使えなかった可能性が。ただ、このようなネットワークのトラブルは今回初めて北大を利用したため起きたことで、ノウハウの蓄積で改善できるはず。といっても次の北大利用はない気が激しくします。あと勉強会の最中私が部屋の外で何かやってましたが、パブリックにできない事なので書きません。

次の北大利用はない気がするのは、道外参加者が多いためです。旅費に比べたら参加費は微々たるものなので、有料の会議室押さえた方が運営側の負担が少ない。

その他色々トラブルもあり、偶然上手くいったこともあるのですが、総じて言えば良く出来ただろうと言いたいです。 もちろんトラブルをゼロにするために、事前の準備をより綿密にする方法はありました(例えば、勉強会当日より以前にゲストアカウントを試しに発行してもらってテストするとか、懇親会の会場は狭くないか下見しておくとか）。ですがやはり有志のボランティアで動かしているので準備に割ける時間は限られているわけです。すでにUstreamの事前テスト(使えなかった。おそらくポート制限かかってる。配信をYoutubeliveにしたのはそのため)とかに運営メンバーの皆さんは相当時間掛けていましたので、トラブルがあっても運営スタッフはよくやってくれたと言わざるを得ないです。皆様お疲れさまです。当日の会場設営に協力して頂いた方々、機材の貸し出しをして頂いた方々もありがとうございます。

勉強会の発表ですが、午前中は私ネットワークのトラブルに対応してたので残念ながらチラチラとしか見てないです(この時まだ申請ミスの連絡を受けてなかった)。昼食時間は、北大の外に出る人もいると考えて時間を多めとっていたのですが、北大の食堂で済ませた人が結構いたようです。午後の発表もトラブル対応してたのであまり見れてないのですが、hgodaiさんのboostのvectorの話とか面白い発表を聞けました。さすがBoost勉強会は濃い話が出ますね。hotwartermorningさんが紹介したbalorは私は使うことはないと思います。。たしかにQtはmocでプリプロセスしてますが、最近のQtはそんなに気にするほどではないような。QtとSTLが相性良くないのはわからなくもない。でもbalorにそれ以上の利点があるとも思えない。江添さんの発表は生で見られて感動。私は数値計算やってるのでアクセスの速いスタックからアロケートできるのは興味深い。ですが私の使わせて頂いてる某計算機はたしか主記憶容量40Tです。スタックってどれだけとれるんでしょうね。懇親会でたしかhgodaiさんから聞いたのですが、スタックからアロケートするとキャッシュヒットしやすいそうです。だからアクセス速いんだとか。後で気がついたのですがキャッシュに乗り切らないサイズだとそれほど意味ない？LTは私頑張って用意できるかなーと思ったのですが妄想だったようです。

懇親会、みんな北海道限定サッポロクラシック飲みましたか？頼まなかった人は何しに札幌来たんでしょう。北海道の地酒もありましたけど、あまり頼んでいなかった模様。懇親会、２次会と盛り上がって、結局カラオケオールでした。Flastさんの近くに行くとカラオケしないで議論が発生するのは前回の教訓として知ってるので近づかない。私は色々やること溜まってて、途中からPC開いて若干作業始めました。最初は見てるだけでしたけど、みんな疲れ始めた頃に歌い出す。そして今期のアニメOPを片っ端からいれる。mikiemon_hさん（だったよね？）が美声で歌ってるのにつられて朝までいてしまった。。また機会があれば行きましょう。

ということで、楽しい勉強会でした。道外から来ていた方々、海を越えて来ていただきありがとうございます。

Matplotlibで軸目盛のフォントサイズを変える

get_xticklabelsを使えばTextオブジェクトを取れる。

あとはTextオブジェクトを自由に設定すればOK。色も変えられる。

import numpy as np
import pylab as pl

pl.figure(figsize=(4,3))
ax = pl.subplot(111)
x = np.linspace(0,2*np.pi)
y = np.sin(x) + -1e-1*x*(x-x[-1])
ax.plot(x,y, lw=2)
pl.xlim(x[0],x[-1])
for i,item in enumerate(ax.get_xticklabels()):
    fontsize = i*2+10
    item.set_fontsize(fontsize)
    item.set_color("red")
pl.show()

ctagsとtaglistをLatexに対応させる

いわゆるExuberant ctagsはコード解析とかに使われる。vimだと関数が定義されている場所にジャンプできたりして便利。vimプラグインのtaglistを使えばctagsで作ったタグ一覧を表示させて、その一覧から選んでジャンプできる。でもLatexではデフォルトで使えない。

設定を追加すれば使える。まずはctags。ホームディレクトリの~/.ctagsに

--langdef=latex
--langmap=latex:.tex
--regex-latex=/\\label\{([^}]*)\}/\1/l,label/
--regex-latex=/\\chapter\{([^}]*)\}/\1/c,chapter/
--regex-latex=/\\section\{([^}]*)\}/\1/s,section/
--regex-latex=/\\subsection\{([^}]*)\}/\1/t,subsection/
--regex-latex=/\\subsubsection\{([^}]*)\}/\1/u,subsubsection/
--regex-latex=/\\chapter\*\{([^}]*)\}/\1/c,chapter/
--regex-latex=/\\section\*\{([^}]*)\}/\1/s,section/
--regex-latex=/\\subsection\*\{([^}]*)\}/\1/t,subsection/
--regex-latex=/\\subsubsection\*\{([^}]*)\}/\1/u,subsubsection/

を書く。次にvimのtaglistプラグインの設定。

let tlist_tex_settings = 'latex;l:labels;c:chapter;s:sections;t:subsections;u:subsubsections'
set iskeyword=@,48-57,_,-,:,192-255

参考 Simple latex ctags and taglist | Caffeinated Code

iskeywordは上のリンクでも軽く言及されてるけどこちらに解説がある

'iskeyword'とは何なのか - LeafCage備忘録

Latex subfigureの番号書式にカッコをつける

subfigureで付けたlabelを本文中で参照しようとすると

図 1.1a

となる。これの正しい呼称はしらない。ラベルフォーマットだろうか？ これを

図 1.1(a)

としたい。プリアンブルにこう書く。

\captionsetup[subfigure]{labelformat=simple}
\renewcommand*{\thesubfigure}{(\alph{subfigure})}

参考: http://www.peteryu.ca/tutorials/publishing/latex_captions

\renewcommandと\renewcommand*との違いが気になった。ググラビティが低くてつらい。

• http://tex.stackexchange.com/questions/1050/whats-the-difference-between-newcommand-and-newcommand
• http://tex.stackexchange.com/questions/86449/what-is-the-difference-between-newenvironment-and-newenvironment

これと同じだろう。