=begin
= 気象解析ライブラリの開発方法
塚原大輔, 堀之内武
最終改訂: 2005 年 03 月 23 日
((:
:))
== 目次
((::))
<<< index.hindex.rd
((:
:))
((:
:))
== はじめに
本チュートリアルでは
GPhys::EP_Flux 内部で用いられている微分演算モジュールを題材に
NArray および GPhys を利用した
Ruby による気象解析モジュールの開発の仕方を
順を追って紹介していきます.
各節の末尾に関連する小技(TIPS)へのリンクを提示します.
気象解析ライブラリを作る上で非常に便利な技などを
紹介していますので併せて参照ください.
TIPS は本チュートリアルの最後尾にまとめて掲載しています.
なお, 紹介しているコードは上記モジュールの一部を
チュートリアル用に簡略化したものです.
== モジュールの仕様を決める
まずモジュールの仕様を決めます.
ここでは以下の仕様のモジュールを作成することにします.
* デカルト座標での微分演算モジュール.
* NArray を対象
* 作成するメソッド
* cderiv(中央差分)
* 引数: データ配列(NArray), 軸配列(NArray), 対象次元(Numeric)
* 返り値: 微分したデータ配列(NArray)
モジュールの空間名は NumRu::Derivative とします. NumRu は
NUMericalRUby の略で, 数値計算や可視化関連のライブラリを
入れるモジュール空間です.
== NArray 版の作成
=== メソッドを定義
# ===== メソッドを定義
#+ メソッドを定義
仕様にしたがってメソッドを定義します. 以下は NumRu::Derivative の 0 版です.
ユーザが呼び出すメソッドは cderiv です.
与えられた配列の両端の境界を 1 グリッド拡張して値を線形補完し(b_expand_linear_ext),
ついでデータ配列, 軸配列の双方の中央差分を計算します(cdiff).
最後にそれらの商を返します.
(()):
((::))
((*TIPS*))
((<例外処理 - 引数のチェック>)) /
((<配列の次元の指定方法>)) /
(())
=== テスト
# ===== テスト
#+ テスト
上記で定義したメソッドをテストします. テストに際して以下の点をチェックすることにします.
* 演算精度
* 解析解との差
* 格子点数依存性
* 不等間隔格子に対する演算結果
* 多次元配列への対応
* 引数チェック時の例外処理
以下にテストスクリプトを示します. 実際の製作過程では一項目ごとにテストを定義, 確認する
方が確実です. まとめて実行するとどこにバグが仕込まれているのかわからな
くなってしまいます.
(()):
((::))
上記スクリプトをコマンドラインから実行するとテスト結果が出力されます.
ただし, モジュールを定義したファイル(derivative_ver0.rb) がカレントディレクトリになくてはなりません.
% ruby test_derivative_ver0.rb
**** equally spaced grid ****
**** single-D ****
dfdx - analytic (except boundary): [mean] 0.03922348996 [max] 0.06451071621
dfdx - analytic (except boundary): [mean] 0.0100170063 [max] 0.01636835692
error change from nx=11->21: 0.2553828409
**** multi-D ****
0.0
test exception successful
**** non-uniform grid ****
"x(11):"
NArray.float(11):
[ 6.28319, 5.68526, 5.14424, 4.6547, 4.21175, 3.81094, 3.44829, 3.12014, ... ]
"x(21):"
NArray.float(21):
[ 6.28319, 5.97675, 5.68526, 5.40799, 5.14424, 4.89335, 4.6547, 4.42769, ... ]
dfdx - analytic (except boundary): [mean] 0.01898250849 [max] 0.03216015526
dfdx - analytic (except boundary): [mean] 0.004934445822 [max] 0.01194461777
error change from nx=11->21: 0.2599469835
上記の結果から, 格子点数が倍になると誤差は約 0.25 倍となり, 定義した微分メソッドは
基本的に 2 次精度であることがわかります. 多次元配列に対する処理も成功しており,
引数のチェックも正しくなされています.
なぜか不等間隔格子配列に対する結果の方が誤差は小さいですが,
たまたま与えた配列が相性が良かったと考えられます.
以上の結果からテストは成功です.
新しいメソッドを付け加える場合はその都度テストを追加していきましょう. また既存のメソッドに手を加えた場合
あわせてテストも書き換えます. テストは頻繁に, かつ慎重に積み重ねることが重要です.
((<モジュールファイルにテストを埋め込む>)) /
((<ユニットテストフレームワーク を用いたテスト>))
=== ドキュメント記述
# ===== ドキュメント記述
#+ ドキュメント記述
定義したメソッドのテストが成功したらドキュメントを書きましょう.
Ruby ではスクリプトファイルの中に埋め込む事を意図して定義
されたドキュメントフォーマット RD が存在します.
今回もスクリプトに RD でドキュメントを埋め込んでしまいます.
RD の書式については((<ドキュメンテーションリンク集|URL:http://www.gfd-dennou.org/arch/ruby/doc-link-j.htm>))
の RD の項のリンクを参照ください.
以下にドキュメントを埋め込んだソースを示します. ここでは英語で書いてますが, 日本語でも構いません.
(Ruby ユーザの間では英語のドキュメントを用意しておくというのが通例のようです.)
(()):
=end
=begin html
=end html
=begin
コツは一つのメソッドが完成するごとに記述していくことです.
ちまちまと少しずつ付け足してます.
メソッドの仕様を変更するたびに更新していくのも重要です.
モジュール全体が完成してから書こうと思うとあとで大変な目に会います.
(更新内容を忘れてしまうという以上にモチベーションが上がりません.)
以上をもって NumRu::Derivative ver 1 の完成です. 今後は新しいメソッドを追加したリ,
既存のメソッドを改良したりしてモジュールのブラッシュアップを図ります.
( 製品版の微分演算モジュールでは
cderiv に加えて 3 点利用二次精度差分(threepoint_O2nd_deriv) が
定義されています. )
その都度, 上記の工程( メソッド定義/改良 --> テスト --> ドキュメント記述) を繰り返します.
== GPhys 版 の作成
以上で, 微分メソッドはできました. しかし, このままでは, 座標とデータは
一まとまりになっておらず, 単位も名前もありません. そこで GPhys 版を作
りましょう. GPhys オブジェクトを対象に微分を行い, 結果もGPhysとして返
します. 微分すると一般に単位も変わりますので, 正しく自動更新される
ようにしましょう. こうすることで, 自己記述的なデータから自己記述的な
計算結果を一行で得ることが出来るようになります. すると, あと 1-2 行で
図にすることも, ファイルに自己記述的な形で出力することもできますので,
とても便利です.
実現方法としては, (1) GPhys にメソッドを追加する, (2) モジュール
関数にする, の2つが考えられますが, ここでは後者を採用します.
前者のほうが自然なのですが, 無闇に自分でメソッド追加をすると,
GPhys 備え付けのメソッドと自前で追加したメソッドの区別が難しく
なりがちという短所もあります. なお, これは GPhys 本家に取り込んで
貰いたいというのが出来ましたら, 積極的に作者まで連絡してください.
(連絡先は((<電脳 Ruby のメーリングリスト|URL:http://dennou-k.gfd-dennou.org/arch/ruby/home-j.htm#mailinglist>))へ)
=== メソッド定義
#+ GPhys の構造の復習
#* GPhys は データ配列(VArray)と格子クラス(Grid) からなる
# * Grid クラスは軸クラス(Axis) を持つ
# * Axis クラスは VArray + 軸情報(配置) を持つ
作成した NumRu::Derivative を利用して GPhys 版 Derivative を作成します.
微分演算のエンジンは NumRu::Derivative にお任せにして,
GPhys 版では物理量の単位と軸名などの属性に留意します.
(()):
=end
=begin html
=end html
=begin
微分演算部分は 22 行目だけです. 特に 27 行目以降はデータの単位と long_name 属性のケア
を行っています.
((* TIPS *))
(())
=== テスト(解析解を用いて)
GPhys 版のテストでは以下の点に注意します.
計算誤差等のチェックは NArray 版の方で行っていますので
ここではしません.
* 属性, 単位
* 微分前後でそれぞれ比較
* 次元の指定
* String/Numeric 双方で指定できるか
以下にテストスクリプトを示します.
長いですが, 解析階も含め自己記述的な GPhys データをゼロから作り出して
るからです. 最初から自己記述的な NetCDF や GrADS ファイル中のデータ
を読んで微分するだけなら, 劇的に短くなります(下記の
((<4.3.3 節|URL:http://www.gfd-dennou.org/arch/ruby/products/gphys/tutorial/ep_flux/#4.3.3>))
参照) .
(()):
((::))
上記スクリプトをコマンドラインから実行するとテスト結果が出力されます.
derivative_ver0.rb, gphys_derivative_ver0.rb がカレントディレクトリになくてはなりません.
最初の行で警告が出ますが, gphys-0.3.5 に NumRu::Derivative および GPhys::Derivative
が既に取り込まれているためです.
% ruby test_gphys_derivative_ver0.rb
./derivative_ver0.rb:8: warning: already initialized constant LINEAR_EXT
******** dimname == 'lat' ********
dfdx - kaiseki_kai (except boundary): 0.01307030336 0.1701432467
**** check attribute ****
name t dt_dlat
long_name temperature d_temperature_d_latitude
units K K.rad-1
******** dim == 0(lon) ********
dfdx - kaiseki_kai (except boundary): 0.01806355861 0.07927459478
**** check attribute ****
name t dt_dlon
long_name temperature d_temperature_d_longitude
units K K.rad-1
******** dim == -1(z) ********
dfdx - kaiseki_kai (except boundary): 3.220621713e-08 7.748603821e-07
**** check attribute ****
name t dt_dz
long_name temperature d_temperature_d_altitude
units K K.m-1
=== テスト(実際のデータを用いて)
ここでは GPhys パッケージに付随するテストデータ(T.jan.nc)を用いて
テストをしましょう. このデータはNCEP 再解析による, 全球の気温の1月の気候値です.
テストデータの詳細は
(())
をご覧ください.
以下テストスクリプトです.
(()):
01 require "numru/gphys"
02 require "gphys_derivative_ver0"
03
04 include NumRu
05
06 t = GPhys::IO.open('T.jan.nc','T')
07 dtdx = GPhys::Derivative.cderiv(t,'lon')
08 newfile = NetCDF.create('dtdx.jan.nc')
09 GPhys::IO.write(newfile, dtdx)
10 newfile.close
T.jan.nc 中の変数 "T" の経度微分を dtdx.jan.nc というファイルに保存しています.
微分がたった 1 行(05 行目)であることに注目です. これだけで
温度場の経度微分が計算できます. 上記のスクリプトを実行して生成された netCDF ファイルの中身を見てみましょう.
% ruby test_gphys_derivative_ver1.rb
% ncdump -h dtdx.jan.nc |less
tcdf dtdx.jan {
dimensions:
lon = 36 ;
lat = 19 ;
level = 9 ;
variables:
float lon(lon) ;
lon:long_name = "Longitude" ;
lon:units = "degrees_east" ;
lon:actual_range = 0.f, 357.5f ;
float lat(lat) ;
lat:long_name = "Latitude" ;
lat:units = "degrees_north" ;
lat:actual_range = 90.f, -90.f ;
float level(level) ;
level:GRIB_id = 100s ;
level:positive = "down" ;
level:long_name = "Level" ;
level:units = "millibar" ;
level:actual_range = 10.f, 1000.f ;
level:GRIB_name = "hPa" ;
double dT_dlon(level, lat, lon) ;
dT_dlon:parent_stat = "Other\n",
"-" ;
dT_dlon:var_desc = "Air Temperature\n",
"A" ;
dT_dlon:missing_value = -9.96921e+36f ;
dT_dlon:long_name = "d_Temperature_d_Longitude" ;
dT_dlon:least_significant_digit = 0s ;
dT_dlon:scale_factor = 1.f ;
dT_dlon:units = "degC.degrees_east-1" ;
dT_dlon:level_desc = "Multiple levels\n",
"F" ;
dT_dlon:actual_range = -72.66724f, -24.35f ;
dT_dlon:dataset = "CDC Derived NCEP Reanalysis Products\n",
"AC" ;
dT_dlon:add_offset = 0.f ;
dT_dlon:precision = 0s ;
dT_dlon:statistic = "Long Term Mean\n",
"C" ;
// global attributes:
:history = "2005-03-11 16:14:03 JST daktu32> NumRu::GPhys::NetCDF_IO.write dT_dlon" ;
}
更に gphys 付属の gpview でコンター図を見てみましょう.
% gpview dtdx.jan.nc@dT_dlon
((::))
このように, 単位などの属性のケアも 1 行でできるようになりました.
=== ドキュメント記述
NArray 版同様にドキュメントをつけて完成です.
(()):
=end
=begin html
=end html
=begin
== 完成
以上で微分演算モジュールが完成しました.
実際に GPhys::EP_Flux の中で用いられている微分演算モジュールは
上記の cderiv に加えて不等間隔格子に対して 2 次精度の差分を与える
メソッド(threepoint_O2_deriv) が加わっていますが, 基本的には上記で示したとおりです.
より複雑なライブラリの作成も基本的にはここで紹介した流れで開発を進めていきます.
意外と簡単だと思われたのではないでしょうか?
自分でもできそうだ, と感じた人はぜひ新たなライブラリを作成してみてください.
繰り返しになりますが, GPhys 本家に取り込んでもらいたいものができましたら
積極的に作者まで連絡してください.
また, 完成していなくてもこんなもの作ってます, という PR も大歓迎です.
((:
:))
== (付録) TIPS
ここでは本文で紹介したライブラリコード中で登場する小技(TIPS)について
解説をします. これらの TIPS は気象解析ライブラリを作成知る上で非常に
便利ですので, ぜひとも参考にしてください.
=== 例外処理 - 引数のチェック
4.2.1 節で提示したコードの各メソッド定義部の先頭では引数のチェックが行われています.
これは想定外の引数が渡された場合は「例外」を上げるというものです.
Ruby では例外処理が言語仕様に組み込まれています.
組み込みメソッド raise は例外を発生させる命令で, 以下の書式にしたがって記述します.
raise 例外クラス, 警告文(String)
ここで例外クラスには引数に関する例外クラス(ArgumentError)などを入れます.
ちなみに例外クラスは省略可能で
raise 警告文(String)
のようにも書けます. この時は RuntimeError が上がります.
とある条件の時に例外を挙げる場合は
if 条件式
raise HogeError, "herohero"
end
や
raise HogeError, "herohero" if 条件式
のように書きます.
より詳しく知りたい方は((<こちら|URL:http://ruby.gfd-dennou.org/tutorial/rakuraku/exception.html>))
をご覧ください.
=== 配列の次元の指定方法
4.2.1 節で提示したコードでは以下のような配列の次元の指定の仕方が多用されています.
これは配列 zの dim 次元目の先頭の値を取り出すものです.
42 val0 = z[*([true]*dim + [0] + [false])] # 配列の先頭の値
これは以下の表記と同じです.
val0 = z[true, tru, true, ..., 0, false]
^^^^^^^^^^^^^^^
true が dim 個
ここで true はその次元の配列要素全てを表しています.
また false は任意個の true を表す「rubber(ゴム)次元」です.
Array どうしの足し算は両辺をくっつけ, 掛け算は個数分の繰り返しであること,
そして最後の引数が * のついた配列である場合, その中身を展開することを利用してます.
( [dennou-ruby:002026] 参照)
ここで Ruby の配列要素の指定の仕方について触れておきます. Ruby の配列は先頭から
0, 1, 2... と数えます. 一方後ろから数える場合は -1, -2, .. となります.
ary = [2, 3, 5, 9, 10]
p ary[0] #=> 2
p ary[-1] #=> 10
p ary[1..-2] #=> [3, 5, 9]
=== NArray 同士の 2 項演算
NArray 同士の 2 項演算では長さ1の次元は任意の長さに自動的に拡張されます.
具体例をお見せしましょう. 以下は [3,3] の配列と [3,1] の配列の和の
演算です. この場合 [3,1] の配列の 2 次元目は自動的に長さ 3 に拡張
されます.
Narray[ [0,1,2],
[3,4,5], + NArray[ [ 0,100,200 ] ]
[6,7,8]]
---> # 内部で自動的に拡張
Narray[ [0,1,2], NArray[ [ 0,100,200 ]
[3,4,5], + [ 0,100,200 ]
[6,7,8]] [ 0,100,200 ] ]
---> # 結果
Narray[ [0,101,202],
[3,104,205],
[6,107,208]]
=== モジュールファイルにテストを埋め込む
4.2.2 節で示したテストスクリプトはモジュールを定義したファイルと独立にしましたが,
短いスクリプトならばモジュール本体のファイルに埋め込む方がよいでしょう.
無駄にファイルを増やさずに済みますし, モジュールの変更に応じてテストを書き
換えたり実行するのが簡単です.
テストを埋め込むにはテストコードを
if $0 == __FILE__
end
で挟みます. $0 は ruby インタプリタによって実行されている
ファイルを表す変数, __FILE__ は自身のスクリプトファイルを表す定数です.
こうするとモジュールファイル自身を実行した場合のみ, テストコードが実行されます.
モジュールファイルを require してもテスト部分は無視されます.
=== ユニットテストフレームワーク を用いたテスト
Ruby にはユニットテストを支援するためのクラス RubyUnit(TEST::Unit) が存在します.
元々外部ライブラリだったのですが, Ruby 1.8 から標準添付ライブラリとして
TEST::Unit というクラスになりました.
個人的な感想として数値計算モジュールのテストにはあまり向かない
(ex. 精度の確認など) と思っているのですが, 興味のある方は利用して見てください.
そして数値計算向けのテストの方法を編み出した方, 積極的に申し出てください.
* (())
* Ruby 1.6 系を利用の方向け
* (())
* (Ruby 1.8 系を利用の方向け)
* (())
* RubyUnit 開発者が書いた解説本. 内容はやや古い.
=== VArray/GPhys の作り方
4.3.2 節で示したコードの 24-25 行目では微分結果の NArray を元に VArray, そして GPhys を構成しています.
24 v_dgpdx = VArray.new(n_dgpdx, v_data, name) # 微分結果の VArray を構成
25 g_dgpdx = GPhys.new( gp.grid_copy, v_dgpdx ) GPhys を構成
VArray は基本的にデータ本体の NArray に, そのデータの名前と属性を加
えた配列です. VArray クラスのインスタンスを生成するには
VArray.new( データ本体(NArray), 属性(Attribute, Hash, or VArray; 省略可, 変数名(String; 省略可) )
のようにします. 上記のように第二引数に VArray を与えると, その属
性がコピーされて, 新たに生成される VArray オブジェクト
に引き継がれます. 属性を陽に定義したいときは Hash を
使うと良いでしょう. なお, NumRu::Attributeはgphysのパッケージ内
で定義されています(解説は略). 第二, 第三引
数ともに省略可能で, 変数名は "noname" となります。nilを与えても
省略したことになります.
上記の例ではデータ本体を微分結果の NArray オブジェクト, 属性は微分前の属性を引き継ぎ新たに
生成した name を用いています.
GPhys はデータ配列(VArray) と格子情報(Grid) からなります.
GPhys クラスのインスタンス生成するには
GPhys.new( 格子情報(Grid), データ配列( VArray or VArrayComposit ) )
のようにします. 先の例(25 行目)では 24 行目で構成した VArray をデータ配列して
元の GPhys オブジェクトのグリッドをコピーして( gp.grid_copy )新たな GPhys
オブジェクトを生成しています.
=end