KuboLog 2005-12-30

2005 年 12 月 30 日 (金)

0740 起床．朝飯．コーヒー．

屋久島常緑樹葉生残過程モデリング，昨日のつづき． Perl による数値計算はおそいので，葉っぱ生き死にの尤度だけを C の関数でやらせよう，という方向にむかいつつある． Perl から C 関数の呼出しに SWIG (simplified wrapper and interface generator) 使ったらいいのでは，と考えて勉強中．

たとえばこの Tutorial にある例にならって calc.c なるものを作ってみる．

double z = 3.0; /* ぐろーばる変数 */
double calc(double x)
{
	int i;
	double y;
	i = 2;
	y = i * z + x; 
	return y;
}

これに対応する SWIG interface である calc.i をこう定義する．

%module calc
extern double z;
extern double calc(double x);

こいつらをコンパイルする Makefile を書いてみる．

TARGET=calc
CORE=/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/
all:
	swig -perl5 -module ${TARGET} ${TARGET}.i
	gcc -fpic -c ${TARGET}.c -I${CORE}
	gcc -fpic -c ${TARGET}_wrap.c -I${CORE}
	gcc -shared ${TARGET}.o ${TARGET}_wrap.o -o ${TARGET}.so

しかし， make するとエラーがでる．

gcc -c calc_wrap.c -I/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/
In file included from /usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/op.h:484,
                 from /usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/perl.h:2344,
                 from calc_wrap.c:291:
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/reentr.h:611: error: field `_crypt_struct' has incomplete type
In file included from /usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/perl.h:3552,
                 from calc_wrap.c:291:
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:199: error: 文法エラー before "off64_t"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:201: error: 文法エラー before "Perl_do_sysseek"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:201: error: 文法エラー before "off64_t"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:201: 警告: data definition has no type or storage class
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:202: error: 文法エラー before "Perl_do_tell"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:202: 警告: data definition has no type or storage class
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:1308: error: 文法エラー before "Perl_PerlIO_tell"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:1308: 警告: data definition has no type or storage class
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:1309: error: 文法エラー before "off64_t"

ということで，計算機の訴えを注意ぶかく検討し，自分の置かれた環境にあわせた対応が必要になる．

いろいろ調べてみてわかったんだけど，いちばん安直な解決策としては，各自の環境において perl -e 'use Config; print $Config{ccflags};' が生成する Perl まわり用コンパイラーオプションをそのまま流用せよ，というもの． Makefile をこう書き直すと無事にコンパイルできる． calc.pm や calc.so が生成される (dir)．

TARGET=calc
CORE=/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/
OPT=-D_REENTRANT -D_GNU_SOURCE -DTHREADS_HAVE_PIDS \
	-fno-strict-aliasing -D_LARGEFILE_SOURCE \
	-D_FILE_OFFSET_BITS=64 -I/usr/include/gdbm
all:
	swig -perl5 -module ${TARGET} ${TARGET}.c
	gcc -fpic -c ${TARGET}.c -I${CORE} ${OPT}
	gcc -fpic -c ${TARGET}_wrap.c -I${CORE} ${OPT}
	gcc -shared ${TARGET}.o ${TARGET}_wrap.o -o ${TARGET}.so

このようにして自作した calc package を呼びだす Perl スクリプト例はかくのごとし．まあ，さすがに簡単さをウリにするだけはある．

use calc;
print &calc::calc(3.0)."\n"; # 9.0 と表示される
$calc::z = 9.0;              # ぐろーばる変数をいじる
print &calc::calc(3.0)."\n"; # 21.0 と表示される

面倒なのは calc::calc() にスカラー以外のオブジェクトを渡したい場合だ．たとえば配列を引数にしたい，という問題を考えよう． SWIG の機能をもう少し使いこなす必要がある．いろいろな解決策がありそうなんだけど，ここではヘルパー関数による方法を試みてみよう． SWIG 本家サイトにある解説文書 Perl Extension Building with SWIG の ``3.4 Helper Functions'' まわりそのままである．

calc.c をこう変えてみる．

double z = 3.0;
double calc(double x, int imax, double* p)
{
	int i;
	double y;
	y = z + x;
	if (imax > 0) {
		for (i = 0; i <= imax; i++) {
			y += p[i];
		}
	}
	return y;
}

対応する SWIG interface calc.i に double 配列の (いわば) constructor / destructor / setter / getter を定義するわけだ．

/* calc.i */
%module calc
%inline %{
double* double_array(int size) {
	return (double *) malloc(sizeof(double) * size);
}
void double_destroy(double *a) {
	free(a);
}
void double_set(double *a, int i, double val) {
	a[i] = val;
}
double double_get(double *a, int i) {
	return a[i];
}
%}
extern double z;
extern double calc(double x, int imax, double* p);

Makefile は上記と同様でよい．こいつの Perl ドライヴァーはこういうかんぢで．

use calc;
sub create_array
{
	my $len = scalar(@_);
	my $ia = calc::double_array($len);
	for my $i (0 .. ($len - 1)) {
		my $val = shift;
		calc::double_set($ia, $i, $val);
	}
	return $ia;
}
$calc::z = 9.0;
my $imax = 3;
my $p = create_array(1, 2, 3, 4);
print calc::calc(3.0, $imax, $p)."\n"; # 22 と表示される
calc::double_destroy($p)

ついでに $p はどういう reference になってるんだ，と表示させてみると _p_double=SCALAR(0x804c464) となってた． C ポインターへの reference?

とりあえず配列を渡すことができれば，葉っぱ生き死に尤度は計算できそうだ． SWIG はなかなかこったつくりになっていて，さらに複雑なデータ構造のやりとりもできるのもウリらしい (けど今回は使わない)．

1040 自宅発．曇． 1055 研究室着．昨晩命じておいた葉っぱ生き死に階層ベイズモデルの MCMC 計算， 10 時間 40 分かかって終了していた (これはまだ all Perl)． 100 MCMC step とばしの 1500 sample の最初の 500 を捨てて作図．収束はすごくおそい． 500 ではだめかも．

このモデルでは休眠などの効果をまだ考慮してないんで，葉齢は「みかけの葉齢」である．で，こいつは「みかけ」であるためなのか，ばらついている，と．

明るさ依存性は「明るいと短命」な方向だけど，よくわからん，と．樹種差・個体差がそろってしまっているのはなんか意味あるんだろな．

(葉面積あたりのとゆー) 窒素割り算値依存性は …… 妥当というかナゾというか．全体としては効果ゼロ，そのばらつきはおそろしく小さい．そして樹種差・個体差が激しい …… ように見える．樹種差・個体差はどちらかといえば，「窒素長命化」な方向で．なンかの交絡なのか?

よくわからん．そもそもこの計算だってアテになるのかアヤしいもんだ．とはいえ改善案もないので，窒素濃度を重量ベイスにして計算やりなおし命じる．無策だ．

SWIG まわりのことなどうだうだ書く．塩寺さんに熱帯林葉っぱデータの準備できたと言われてあせる．しばらくお待ちください，とゆーことで．

昼飯．今日は院生密度がさすがに低い．

1/2 にいきなり京都往復，というよーなハナシが．うーむ．とりあえずいろいろ調べてみる．

のんびりと修論にとりくむ石橋君のミズキモデリングこんさる．個体差が見えにくい．

歳末も遺伝子をよむ亀山さんとツガザクラ雑談．この植物，クローナル成長が重要なのか重要でないのか，まったくわからなくなった．場所とりは何で決まるか …… さらにこれと自殖の可否との関係もよくわからん．たとえば，単一クローンにみえるのが全部自殖種子由来，とかありえんのかしらん (親がヘテロなんだけどどちらの染色体も劣性致死になってる，とか)．自殖率と近交係数から近交弱勢が決まる，と教えてもらったんだけど，そうなると近交係数の空間構造依存性が気になってきたり …… このあたり，あまり重要とは考えられてないよーな．

SWIG による屋久島モデリングの尤度計算高速化，つまり Perl を部分的に C に置き換える，という試み．やってみたら， 100 MCMC step 計算に要する時間が 56.2 秒から 36.3 秒に短縮された …… うーむ，改善と言えるかどうかびみょーな気も．もうちょい速くなってほしい．

尤度計算だけでなくメトロポリス法とかでもてこずっているのかね …… よくわからんな．ぎぶすさんぷらーだと速くなる，とか? どうなんだろう．

ぎぶすさんぷらー化するのはたいへんなので，安易なトコロから …… 正規乱数生成とか事前分布まわりの尤度計算も C に置き換えてみる． 100 MCMC step 計算に要する時間が 36.3 秒から 35.3 秒に短縮された …… うーむ，努力がむくわれん． Perl もそこそこに速い，ということかな (いちおーコンパイルしてるそーで)．

まあ，これ以上速くしようとするなら，コード全部を呪われ言語 C++ で書き直すしかあるまい．その時間的よゆーがなく，無理っぽいような気がする．

ぢたばたしても改善されないので，撤退準備 …… すると，なぜか晦日も R 作図にはげむ西村さんから「par(mfrow = c(縦, 横)) で図を並べるとつぶれてしまう」質問．これは FAQ というか …… (grid を用いない場合に) 各わく内でつぶさないポイントは par(mar = c(下, 左, 上, 右), mgp = c(軸名, 軸数字, 軸線)) の調節 (特に 上, 右)，そして出力 device のサイズ調節である．

2210 研究室発．買いもの．食料の多くを輸入にたよっている北海道では，冬の荒天による輸送とどこおりに敏感である …… 野菜の値段とか，が．たとえば，ホウレンソウ一把 238 円．私は 1-2 食ぶん 150 円を超える野菜はなかなか買えない． 2230 帰宅．晩飯．

窒素の重量ぱーせんととやらを使ったモデルの計算結果がでていた． 501 から 1500 番のサンプルを使って様子をみてみると，昨日朝に「とりあえず」作図したのとあまり変わらないような気がする．

おそらく窒素まわりが明るさと交絡してるンでは，という予感．シュート伸長の休眠とかをいれると，何か変わるんではないかな．現時点の予想としては，「明るさ」の効果が明確になり (暗いところでよく休眠するので)，窒素うんぬんはそれによって抹殺される …… となってればハナシは簡単なんだけどな．
とはいえ，このあたりの設計・実装にはもうちょい時間かかりそう．明日ぐらいには何とかなるか? ともあれ，現在うごかしている MCMC 計算プログラムの信頼性には，まだまだ疑うべき点が残されているように思うので，昨晩から本日にかけての計算の検算 (異なる初期状態からの MCMC 計算) を A 棟 8F の Dell 機に命じておく．
[今日の運動]
- うーむ，このうんどう不足状態はどうやったら解消されるのか …… ぢみちに腕立伏せとか? やらないよりはよほどマシだろうな (肩こりにも効きそうだし) ……
[今日の食卓]
- 朝 (0820): 米麦 0.7 合．ダイコン・ニンジン・ネギ・ナメコ・タラすりみ味噌汁．
- 昼 (1330): 研究室お茶部屋．こんびに握り飯二個．
- 晩 (2250): 米麦 0.7 合．ダイコン・ニンジン・ネギ・ナメコ・タラすりみ味噌汁．

本日 (kubolog20051230) | 次の日 | 1 日前 | 7 日前 | 31 日前 | 365 日前 | top

KuboLog | KuboWeb

ぎょーむ日誌 2005-12-30

2005 年 12 月 30 日 (金)