KuboLog 2005-08-17

2005 年 08 月 17 日 (水)

0820 起床．朝飯．コーヒー．洗濯． 0935 自宅発．晴． 0950 研究室着．

昨日からてぬき光計算の具体的な手順について検討している．一週間前に示したように，基本アイデアはこういうかんぢで．

方向と重みの関係は一昨日のこれ．

[方向と重み]

光量のもっとも単純なモデルを考えた場合，面積あたりのエネルギーは sinφ に比例する (φ は「視線」の仰角)．かかる格子状配列の世界だと

sinφ = 
		sqrt(dz² /
		(dx² + dy² + dz²))

と簡単化される (d* は視点座標と視線上点の座標の差)．

昨晩の北大構内走中，あるいは今朝おきぬけのぼーっとした状態で列挙した手順の断片あれこれ．

明るさ減衰計算は全部 log(明るさ) でやる --- これによって正しく相乗平均的に計算できるだけでなく， Fbox ごとに計算をうまく分割できる
各 Fbox が foliage vertical distribution (FVD) をもつ
FVD への加算・減算関数をつくる --- とりあえず枚数の増減だけがわかればよいはず --- いや，このあたり再検討が必要
FVD から三通りの対数減衰速度を計算しておく --- 加重ナシ (自分用)，上加重 (他人用)，下加重 (自分用)

…… ありゃ．意外とこれだけで MCMC 駆動に必要なてぬき明るさ計算は完結してるのかしらん?

さて，こういう計算専用のクラスを作るかどうか，だが．うーん …… とりあえず LightProcessor クラス，とゆーのを作ってみるか．

まずてはじめに， Fbox に全部おしつけてた作業の一部を LightProcessor に分担させればいいだけなんだが …… なぜかばぐとりにちょっと手間どる．仕事のひきつぎを確認．時刻はすでに 1125．

次第にタテわり官僚主義になりつつある母子里 Gibbs 林冠計算プログラムの構造改革のつづき．じつは光計算にまつわる仕事は， LightProcessor の「上司」である Fbox をトバしてもぎょーむがススむ場合が多い (つまり下っぱの LightProcessor どうしで情報交換すればよい)，とわかってきたんで，そのように「風とーしのよい」しくみに変更していく．変更 & 試験運転のくりかえし．

てなことやってるうちに，上の図のような局所明るさ計算に入りこむ． LightProcessor で処理すべき計算はぜんぶ書けた．時刻は 1355．昼飯にするか．

北大生協に昼飯調達にでる．研究室にもどって昼飯．

またプログラミングのつづき．昼飯前の作業がうまくできてるのを確認する試験運転．次にこの計算世界でいちばんエラい Canopy から LightProcessor を順番によびだして明るさ計算やらせるところを作る．

ケガれ言語 Perl のおぶぢぇくと指向なデストラクターって使ったことなかったんだけど，今回は使ったほうがよいかも，ということでデストラクター試験コードを作ってみる．

#!/usr/bin/perl -w
#-------------------------------------------------------------
package Leaf;
use strict;
sub new # コンストラクター
{
	my ($class, $id) = @_;
	my $self = { 'id' => $id };
	bless $self, $class;
	printf "Leaf %i is constructed.\n", $self->{'id'};
	return $self
}
sub DESTROY # こう書くとデストラクター関数として処理される
{
	my ($self) = @_;
	printf "Leaf %i is destructed.\n\n", $self->{'id'};
}
#-------------------------------------------------------------
package main;
use strict;
my %hash_leaf;
$hash_leaf{'1'} = Leaf->new(1); # hash に放りこんで……
delete $hash_leaf{'1'}; # 削除 
my @array_leaf;
push @array_leaf, Leaf->new(2); # array に放りこんで……
shift @array_leaf; # 削除
my $leaf = Leaf->new(3); # プログラム終了時にきえる $leaf
__END__

実行結果はこういうあたりまえなもので．

Leaf 1 is constructed.
Leaf 1 is destructed.
Leaf 2 is constructed.
Leaf 2 is destructed.
Leaf 3 is constructed.
Leaf 3 is destructed.

これを今つくってるほうのプログラムの Leaf でも追加する．生成するときにコンストラクターで LightProcessor に ->Add_leaf($self) して，捨てるときに同じく LightProcessor に ->Delete_leaf($self) する，と．

これで葉っぱの増減にともなう処理を簡単化できるんだけど …… いやはや，私が上のようにしといたことをうっかり忘れて，ヘンなコードを書いてしまいそうでおそろしい．

といった失敗しないように注意コメントを Delete_leaf 関数につけてるうちに，さらに強まった my マヌケさ (use strict に表現してみました) を発見してしまった．ここで Delete_leaf を呼び出せる状況においてはデストラクターは発動しない．なぜかと言えば， LightProcessor 内でその Leaf への参照が保持されているからだ (Delete_leaf はその参照を削除する関数)．

ついでながら …… Perl の garbage collection は mark & sweep 方式．

だめだ． DESTROY の使いどころない．いやはや．やはり葉っぱ捨て処理のときにまず

$light_processor->Delete_leaf($leaf) してから
%{$canopy->{'foliage'}} (蛇足ながらこういうのが Perl のケガれ書法だっ) に格納されている最後の reference を消す

と．ややこしくて，自分でもド忘れしそう …… なんで，ぎょーむ日誌に詳細を記してみる．

とりあえず葉っぱ増減にともなう問題はあとまわしにして，局所明るさ計算の試験運転とばぐとり．

局所明るさに関する試験運転，とりあえずできた．よい結果 & よくわからん結果でる．よい結果は計算時間が短かった，ということ．葉数 10000 枚ぐらいだと最悪でも 1 秒ぐらいで計算がかたづく．これはこの先にひかえるめんどう MCMC 計算をススめるのに助かる．

よくわからん結果はこれ．

昨日の方法でもって「てきとーな初期状態」 (上の図の左) を作って，とりあえずその局所明るさ評価をやってみたんだが，それが右の図．いくら「てきとー初期状態」でもこれは明るすぎるよな．葉数 10000 枚とゆーか LAI が 4 ってのはそんなに小さい数ではないはずなんだが．

計算まちがいか，あるいは手ぬき計算方法が手ぬきしすぎている可能性がある．ちょっと調べなおすか．

えーと，一個のハコ (60cm × 30cm × 60cm) に直径 10cm の球状の葉っぱ押しこんでいく場合 …… LAI = 4 を実現するには 91 個ほどつめこむ必要あり．で，計算してみるとさすがにそれだけつめこめば明るさはゼロちかくまで減衰する，とわかった．まず根本的なところはまちがってない，と確認．

なーんか単純なミス，という気もするんだが ……

しかし本日もケガれ言語プログラミングでアタマが消耗してんで撤退． 1920 研究室発． 1935 帰宅．

しかしやはりヘンてこ計算結果は気になるんで気合いっぱつバグ狩りの解析 …… 30 分で原因特定できた．計算まちがいではなくプログラムのバグだった．

いやはや，上のほうで「風とーしのよい」しくみうんぬんといった仕事分割やった，と書いたけどそのときに入ったエラーだった． Fbox 間の接続がヘンになっていた．かくのごとく，構造改革は難しい．

で，これで明るさ計算やらせてみると …… おー，われながら「うまい初期状態」生成させる方法を発明しちゃったよーで．この段階ですでにかなり「あってる」よーに見える．

えーと上の図はどちらも横軸がダケカンバ林分内の局所明るさの観測値で，縦軸がモデルの「初期状態」葉群配置の仮想林分で実現している局所明るさ (左は常数軸，右は両対数軸)．明るさは 0-1000 であらわされている．対数軸でみると観測値 50 (相対明るさ 5%) より下の観測点はほとんどないんだけど，モデル初期状態ではさらに暗い場所が林内にたくさん存在することになっている．

まあ，たぶん「LAI = 4 ぐらい? だったら葉数は 10000 ぐらい?」などなどといった前提のどこかがまちがっていたんだろう．そしてこの葉群「初期配置」の方法には何の正当化もない．なぜかとゆーと …… 私の直観だけで導出された計算方法であるからだ．科学というよりオカルトですな．実際のところまじめに検討すればいろいろな欠陥を見つけられる．

そこで …… こういった不適切な前提だとか，あるいは一見もっともらしいけどじつは根拠ふりーで決めてしまっている葉の位置なんかを，統計学・統計物理学の基本わざの連鎖わざでもって改善していくのが Markov Chain Monte Carlo (MCMC) 計算なのである．とゆーか，そうであってほしい，ということで今回の問題に取り組んでいるわけだ．

たとえば，葉数に関しては平均 10000 のポアソン分布を (ベイズ統計学ふうにいうなら) 事前分布として与えたとしよう．そのような状況での MCMC サンプリングの結果として生成される事後分布はたとえばもっと低い平均値をもつ確率分布になる，といった「改良」が得られるかもしれない …… そういう研究をやってるわけで．

いやー，「発注者」小林さんが明日から札幌に来られるんだけど (その後，このモデリング対象である母子里の森林に調査へ)，その前に MCMC 計算以外の難所はなんとかすべて片づいた，か．めでたしめでたし．

一件落着したんで 2030 自宅発北大構内走．かなりよろよろ． 2125 帰宅．ばてばて．体重 72.8kg．晩飯．

本日の Gmail アカウント出荷は 1 個．まだまだありますよー

[今日の運動]

北大構内走 55 分間．

[今日の食卓]

朝 (0840): 米麦 0.7 合．キャベツ・モヤシ・ネギ・ブナシメジ・サケあらの味噌汁．
昼 (1440): 研究室お茶部屋．北大生協ツナサンド．
晩 (2230): 米麦 0.8 合．メカブ．ニラ・エノキダケ・卵の炒めもの．キャベツ・モヤシ・ネギ・ブナシメジ・サケあらの味噌汁．

ぎょーむ日誌 2005-08-17

2005 年 08 月 17 日 (水)