2009年6月28日日曜日

R で描画しよう

MathematicaがUpdateしてから調子悪い。
描画するとラベルの字がつぶれるつぶれる。25万もするくせになにこれ、うぇーん。
EPSで出力するとフォントつぶれは回避できるんですが、MathematicaのEPSはクソ重くて、それだけで10Mとか論文に貼り付けるにはあまりに激しい容量です…。
そこで、3Dプロットを美しく書けるフリーソフトを探していたわけです。
要求することはただ一つ。離散的な3次元のデータをデータ間の補間を行いつつ3Dプロットすること。
まぁGnuplotでもそれなりのものはかけて、
set pm3d
でpm3dを読み込んでやれば、splotで
http://t16web.lanl.gov/Kawano/gnuplot/plotpm3d.html
ここにあるような3Dのカラーグラフはかけます。しかし、いかんせんGnuplotは3Dデータの補間がないようなのですよ。csplineは基本2次元plot用だし…。
そこで手をだしてみたのが"R"です。
"R"の存在は前々から知っていてその異様にシンプルな名前にびびって手を出せていなかった代物です。統計解析ツールというかそれ専用の言語+統合開発環境です。言語なのでガチでプログラミングぽいです(ToT)
さて、Rの強みの一つはグラフィックが美しいことです。そして統計解析専用に作られているんで、補間用のパッケージとかもごろごろ転がっています。といいことづくめなようなのですが、めっちゃ躓きました。
グラフィックを出力したいのに、デフォルトのグラフィック用の関数が、等間隔のグリッドにしか対応してねぇ…。
私の場合、計算のデータポイントは基本、変動が大きそーなとこは細かく、どーでもよさそうなところは粗く、対数的にとっているのでこれはかなり困りました。
4時間ほど苦闘した結果、
akima(Akima先生?とにかくありがてぇ)っていうライブラリ中にあるinterpという関数を使って補間をつかって等間隔グリッドを作成、データをリストに格納→描画
っていう流れでなんとかなりそうな予感です。
備忘録として、
CSVからデータを取り込み→それをリストに格納→interp→描画
のテストコードをあげておきます。


book1<-read.csv("book1.csv")
library(akima)
rlist<-as.list(NULL)
rlist$x <- c(book1[,1])
rlist$y <- c(book1[,2])
rlist$z <- c(book1[,3])
rlist.li<- interp(rlist$x,rlist$y,rlist$z,
xo=seq(min(rlist$x),max(rlist$x),length=50),
yo=seq(min(rlist$y),max(rlist$y),length=50))
image(rlist.li)
persp(rlist.li)



read.csv():CSVから読み込み
c():リストへの要素の追加や結合
persp():俯瞰図の描画
interp:補間。lengthでメッシュの数を決めることができます

CSVファイルとして

x y z
-5 -5 50
-4 -5 41
-3 -5 34
-2 -5 29
-1 -5 26
0 -5 25
1 -5 26
2 -5 29
3 -5 34
4 -5 41
5 -5 50
-5 -4 41
-4 -4 32
-3 -4 25
-2 -4 20
-1 -4 17
0 -4 16
1 -4 17
2 -4 20
3 -4 25
4 -4 32
5 -4 41
-5 -3 34
-4 -3 25
-3 -3 18
-2 -3 13
-1 -3 10
0 -3 9
1 -3 10
2 -3 13
3 -3 18
4 -3 25
5 -3 34
-5 -2 29
-4 -2 20
-3 -2 13
-2 -2 8
-1 -2 5
0 -2 4
1 -2 5
2 -2 8
3 -2 13
4 -2 20
5 -2 29
-5 -1 26
-4 -1 17
-3 -1 10
-2 -1 5
-1 -1 2
0 -1 1
1 -1 2
2 -1 5
3 -1 10
4 -1 17
5 -1 26
-5 0 25
-4 0 16
-3 0 9
-2 0 4
-1 0 1
0 0 0
1 0 1
2 0 4
3 0 9
4 0 16
5 0 25
-5 1 26
-4 1 17
-3 1 10
-2 1 5
-1 1 2
0 1 1
1 1 2
2 1 5
3 1 10
4 1 17
5 1 26
-5 2 29
-4 2 20
-3 2 13
-2 2 8
-1 2 5
0 2 4
1 2 5
2 2 8
3 2 13
4 2 20
5 2 29
-5 3 34
-4 3 25
-3 3 18
-2 3 13
-1 3 10
0 3 9
1 3 10
2 3 13
3 3 18
4 3 25
5 3 34
-5 4 41
-4 4 32
-3 4 25
-2 4 20
-1 4 17
0 4 16
1 4 17
2 4 20
3 4 25
4 4 32
5 4 41
-5 5 50
-4 5 41
-3 5 34
-2 5 29
-1 5 26
0 5 25
1 5 26
2 5 29
3 5 34
4 5 41
5 5 50

こんなのを用意しとくと


こんなかんじのができます。
これはx^2+y^2を描画してみたものです。

2009年6月27日土曜日

旧友とランチ@牡丹園

久しぶりに高校時代の友人に会った。なつかしい。メールのやり取りはちらほらあったのだが、会うのはずいぶん久しぶりである。
関東からほぼ日帰り、帰りは夜行バスらしい。母は強しだな~。

昼ごはんを一緒に食べようということで、EST内にある神戸別館牡丹園に行った。本店は名前の通り神戸にあって、昔は千里中央の阪急にも支店が入っていた。ここの焼きそばやひき肉を甘辛く炒めた物をレタスで包む料理が好物で、何度か頼んで祖父母に連れて行ってもらった記憶がある。牡蠣のお好み焼きはとくにハイパーおいしかった記憶がある。十数年ほど前?に千里中央店がなくなってしまい、それ以降神戸はさすがに遠くてなかなかいけなかったのだが、最近ESTに入ったときき、行く機会をうかがっていたのだ。梅田飲食店激戦区なだけあって、ランチメニューやらオンラインクーポンなどあってお得である。

いろいろ話を聞いていると、女の集団のなかに入るってのは大変そうだなぁと思った。女性がある程度多いとなぜか「女王蜂」が誕生してそれ以外が(男女を問わず)「働き蜂」的ポジションでへいこらして持ち上げるてな構図になりがちっぽいなぁ。それが集団を平和に維持するための最適なモデルなのかもしれんがな、個人的にはちょっとなぁ。ネットワーク理論でこういうのも解けるのかしら。

あと、母親に関する思い出で共通する部分があってものすごくウケた。あぁみんなやっぱりそんな体験があるんだなと。

2009年6月15日月曜日

Tex用エディタ

texのエディタでかなり高機能なものを見つけた。
ここ
http://www.juen.ac.jp/math/nakagawa/texguide.html
にある、EasyTex。数学の先生が作られただけあって、数式関係の機能が凄い。
数式の行数を聞いてくれて、左寄せするテンプレートまで作ってくれるとは。腐るほど数式かかなあかん理論系にはありがたや。

Winshellも高機能なんだけど、Vistaだとどうもうまくいかないみたいで、後輩が困ってた。こればVista対応とあるし、勧めてみよう。サクラエディタにマクロで十分便利なんだが、ギリシア記号の入力支援とかラベル参照のときの支援があると断然楽だと思う。

もうひとつ、EclipseのTex用プラグインがなかなか秀逸そうなので、そっちも試そうと思う(タブで開けないと気持ち悪い病にすっかり侵されている)

2009年6月7日日曜日

f77とf90の混在

さて、最近良くも悪くも忙しい。論文なぞ書いていて、ふと、「この数値計算ほんとに大丈夫なんだろうか」という不安が頭をよぎりだしたりして、D論を終わらせる前に計算の妥当性をちゃんと調べておこうという気になった。まーすでに何本かジャーナルに出てるんで妥当でないと困るんだが。

なんでそんなことが気になるかというと、NRGにおけるちょとしたテクがいまの系ではうまく働かないからだ。NRG(数値くりこみ群)は伝導電子バンドを対数的に離散化することで、問題を有限サイズ系の対角化にはめこんでいる。離散化の影響を最低限に抑えるために、いろんなメッシュとろうぜってのを試すとちょっとアレ?ってな感じだったりする。いまの自分のコードでは一般的な場合は影響がほとんどないんだが、伝導バンドとの混成にエネルギー依存をつけると、あやしい。一般的な場合に影響がでてないので、NRGの本体ではなくそこにほりこむパラメータ類の計算が怪しい。

パラメータの類はMathematicaで計算していた。Mathematicaは高性能かもしれんが、ブラックボックスなので不安である。なんでMathematicaかというと、計算に特殊関数のBessel関数とその数値積分がいるからだ。ニューメリカルレシピは単精度だったのでつかえねぇ…。というわけでMathematicaつかってたんだが、つい最近、高精度のBessel関数のフリーなライブラリを見つけた。
http://www.kurims.kyoto-u.ac.jp/~ooura/index-j.html
オオウラ先生ありがとぅぅぅ。Bessel関数だけじゃなく、高精度数値積分のパッケージまで作っておられる。まさにネ申。というかここのパッケージ使えば大半の物理計算ができてしまうと思われ。

恥ずかしながら、DE(二重指数関数型)公式という数値積分法は知らなかった。もう一個のIMT公式もだけど、日本人が発見した方法らしい。日本ってやっぱ数学つよいんだな。変数変換をもちいて区分求積することで猛烈に精度がよくなるらしい。なんでうちの専攻の授業ではガウス積分とかばっかでもっとつかえそうなこっち教えないんだ?

さて、ありがたくもらってきたんだが、ソースはf77で書かれている。私の普段使ってるのはf90。でも大丈夫。基本的にはおなじコンパイラ(g95,ifort)でコンパイルするので、-cオプションでオブジェクトファイルにしておいてリンクすれば.f90のファイルから呼び出すことができる。

ただし、問題が一つ。f90は結構型宣言について厳しい。私は、先頭にimplicit noneつけないと気持ち悪いんだが、これをするとfunctionを使うときにやっかい。たとえば、外部のファイルにfunction hogeがあるとして


program main

implicit none

integer :: i
double precision :: x,res
external :: hoge

res=hoge(x)

end program main

は、hogeが型宣言されてないとしてはじかれる。関数を double precision :: hoge みたいに型宣言してもあいてがf77なので仕様が違う→結構ヘンなことになりそうで怖い。

program main

implicit none

integer :: i
double precision :: x,res

interface
double precision function hoge(x)
double precision , intent (in) :: x
end function
end interface


res=hoge(x)

end program main

とinterface文でhogeの属性を指定してやるのがおそらくベスト。うぅ。いちいち使うたびにinterfaceかくのめんどくさいので、


subroutine ldbesj0(x,res)
implicit none

double precision :: x,res
interface
double precision function dbesj0(x)
double precision , intent (in) :: x
end function
end interface

res=dbesj0(x)

end subroutine ldbesj0

というのを各ベッセル関数について作った。こうすると、ライブラリにしたときに、ldbesj0とかを呼べばinterfaceについて考えなくてもよいしなぁ。というわけで、明日はライブラリ作ろうと思います。