Ruby Library Report 【第 5 回】 数値計算と可視化

著者: 堀之内 武、編集: 立石 孝彰

はじめに

Ruby は数値計算には向かないと言われる。確かに Ruby は速くはない。しかし、 Ruby によりプログラミングが簡単になることは大きな魅力である。数値計算 をする場合、一般には何か入力があって出力がある。両者が少数の数値のみか らなる場合はあまり気を使う必要はないが、ひとたび教科書的な演習問題を離 れ、現実的な問題に乗り出すや否や、入力も出力も単純/少量でなくなり、ま た、似たようなことをするのに複数の形式のデータを扱わねばならぬといったこ ともあろう。数値とともに名前やら日付やらといった文字情報を扱う場合もあ る。このような場合に、洗練された使いやすいオブジェクト指向言語である Ruby は大きな助けとなる。また、数値計算本体においても、することが複雑 である場合、クラスを導入することでプログラムの見通しが良くなり、可読性・ 保守性が上がるケースも多い。つまり、Ruby を数値計算に用いるのは伊達や 粋狂でなく、現実的な問題により良く対処するためなのである。

では、速さについてはどうなのか。Ruby の 1 ステップは、大ざっぱには C の 100 ステップ程度に相当する。しかし、以下のようにそれでも問題が少ない場 合がある。そして、実はそういう場合はかなり多いのである。

  1. 計算量自体が少ない場合。例えば C で 1 秒以内に済む計算が、Ruby では 1 分かかるかもしれない。しかし、そのためのプログラム開発に C では 5 時間かかり、Rubyでは 30 分程度で済むならどうだろう。もしもそのプログラムを実行するのが高々数10回程度としたら、Ruby を使うことはおかしな選択ではない。
  2. 拡張ライブラリーが利用できる場合。特に多量のデータを配列に入れて計算に用いる場合、後述する NArray という拡張ライブラリーを用いることで、ループが C で回るようにすることができる。すると、計算回数に比べて Ruby のソースコードレベルでの処理数は格段に小さくなり、高速に計算を行い得る。

2. の場合も、Fortran や C のみで書く場合に比べればやはり遅いであろう。 しかし、要は数値計算と言えども開発にかかる時間やプログラムの保守性を考 慮すべきであるということである。すると、Ruby の利用は十分現実的な選択 肢の一つであろう。

Ruby は C で実装されており、C で容易に拡張ライブラリーを書くことができ るのは、大きなメリットである。Ruby 組み込みのクラスの多くは C で書かれ ており、それを Ruby で利用できるようにする仕組みがそのまま使えるので、 容易であるだけでなく、移植性も高い。C だけでなく、これまでに多くの数値 計算ライブラリーが開発されてきた言語 FORTRAN77 のライブラリーを利用す ることも、比較的容易である。C や FORTRAN の既存のライブラリーを Ruby に組 み込むことで、過去の資産をより使いやすい形で活かすことができる。うまく いけば、実行速度をほとんど落とすことなく、開発速度を飛躍的に高めること ができる。なお、後述するように、数値計算結果の可視化においても、拡張ラ イブラリーは大きな役割を果たす。

以下ではまず、数値配列のデファクトスタンダードである拡張ライブラリー NArray を紹介する。均一な型のデータをメモリー上の連続領域に配する NArray は、同様のデータ構造をとる多くの C や FORTRAN のライブラリーと Ruby の橋渡しもする。NArray に次いで、可視化ライブラリー RubyDCL を紹 介する。これは DCL という可視化 ライブラリーを「ラップ」する拡張ライブラリーである。 さらに、3次元可視化ライブラリー VTK の拡張ライブラリーもごく簡単に紹介する。 次に、数値計算のための拡張ライブラリーを取り上げる。 代表として、フリーで包括的な数値計算ライブラリー GNU Scientific Library (GSL) を用いる Ruby/GSL を紹介し、さらに歴史のある商用数値計算 ライブラリー SSL2 と Ruby をつなぐ Ruby-SSL2、そして有名な行列計算ライ ブラリー LAPACK ベースの Linalg につ いて触れる。数値計算の入出力に関しては、バイナリー数値データを構造化し て格納する NetCDF の拡張ライブラリーを紹介する。 最後に、応用例として GPhys という格子点データ取り扱いライブラリーを紹介する。 そこで用いられる、物理量の単位取り扱いライブラリーにも触れる。

筆者の専門は大気物理学である。前述の DCL も気象学者が開発したものであ るため、一般的な数値データの可視化機能に加え、データの欠損がありえる物 理的なデータに対応でき、さらに地図投影など地球科学固有の問題に対処する 機能がある。NetCDF も、気象・海洋学でよく使われるデータの形式兼入出力 ライブラリーであり、任意のバイナリー数値データが格納できるのに加えて、 物理量データの構造化ができるのが売りである。しかし、他の分野でも使われ ており、筆者の知る例では Intel に対抗する某半導体メーカーで使われてい る。最後に紹介する GPhys は、それまで紹介してきたライブラリーの上に構 築されたライブラリーである。物理データに関する高度な抽象化を行っており、 筆者が日々大気データを解析する際に直接使うライブラリーがこれである。基 礎となる個々のライブラリーに比べれば、恐らく利用者は限定されるが、本稿 で紹介するライブラリーの活用例として見て頂きたい。オブジェクト指向のう ま味を活かしてライブラリーを作り込むことで、データ処理が大幅に省力化さ れる例となっている。

本稿で取り上げるライブラリーは、「科学」における Ruby の利用という観点 からは、限られたものとなっている。興味を持った読者は、NArray の作者の 田中さんがまとめられている Wiki ページ 「Ruby for Scientific Computing」 を参照されたい。簡単な説明を付した良いリンク集である。

始める前に

本稿で紹介したライブラリーの幾つかは、筆者が関わる 「(地球流体)電脳Rubyプロジェクト」 で管理している。 ここでは、それ以外のライブラリーについても情報をまとめ、さらに Linux、 FreeBSD や Windows 用のソース及びバイナリーパッケージを提供している。 例えば、Windows 用には Visual C++ でコンパイルされた個々のライブラリー に加え、Ruby 本体も含めすべてを一発でインストールできる「setup一発」イ ンストーラがダウンロードできる。Debian 及び Vine Linux では apt 一発で インストールでき、FreeBSD でも"ports" が用意されている。 なお、地球流体電脳クラブ は、大気・海洋など、地球・惑星の流体の研究における計算機利用の諸側面を 扱う有志の集まりである。

レポート NArray

登録データ

概要

NArray は多次元数値配列のクラスである。データ型はオブジェクトごとに均 一で、例えば倍精度実数をもつ NArray オブジェクトでは、すべての要素が倍 精度実数である。FORTRAN など多くの言語の多次元配列と同様、各次元がそれ ぞれ均一な長さを持つ、言わば長方形的な配列である。NArray は、Cで書いた 拡張ライブラリーとして実現されており、データを入れるメモリー上の連続領域 へのポインターと配列の次元数 (ランク) や各次元のサイズなどからなる構造体 を Ruby オブジェクト化するものである。サポートされる型は1,2,4バイト整 数 (1byte の場合のみ符号なし整数)、単精度・倍精度の実数と複素数、そして Ruby オブジェクトである。最後の点は少し補足が必要でろう。実は、NArray は、Ruby の実装における C レベルでの VALUE 型を利用し、任意の Ruby オ ブジェクトを持つ多次元配列を作ることができるのである。ただし、この場合 は C を使うことによる高速化は期待できない。以下では本題の数値型の場合 に限って説明する。

数値型の NArray では、各要素は Ruby オブジェクトにはならず、C ネーティ ブな型で表される。ループも C レベルで回るため Ruby 組み込みの Array に 比べて遥かに高速な演算が出来る。四則演算等は要素ごとの演算として定義さ れている。

 a * b

は [ a[0]*b[0], a[1]*b[1], .. ] を意味するといった具合である。 組み込みの Array 同様、各次元の添え字はゼロから始まり、負の添え 字は後ろから数えることを意味する。

NArray には、NMatrix, NVector という行列・ベクトルに特化したサブクラス がある。後述するように演算が行列・ベクトル用に再定義されており、さらに LU分解等のメソッドがサポートされている。

NArray の次元の並びは FORTRAN 式である。即ち、メモリー上連続なのは、最 初(左端)の次元であり、最後 (右端) の次元がもっともゆっくりとまわる、 column(列)-メジャーな配列である。しかし、実は NMatrix は、第 1 添字が列 で、 第2添字が行と、行列ならぬ「列行」になっている。よって、本来の行 列の意味に立ち返ると、実は列に関して連続で行に関してより「大きく」回る、 row (行)-major であることは注意を要する。このため、添字順に「行列」に する慣例に従っている外部ライブラリーに対しては、むしろ C との相性がよい。

サンプル

初歩

最初に簡単な例を示す。

na1.rb

   1|require "narray"
   2|na = NArray.float(3,2).indgen!         # indgen! --> 0,1,2,..
   3|nb = NArray.to_na( [ [1.0],            # Will be a float NArray
   4|                     [10 ] ] )
   5|nc = NArray.sfloat(12).random!
   6|print "na = "; p na
   7|print "nb = "; p nb
   8|print "nc = "; p nc
   9|print "na + nb = "; p na + nb
  10|print "na >= nb = "; p na >= nb
  11|na[na >= nb] = 999
  12|print "na[na>=nb]=999 --> na = "; p na
  13|
  14|print "na[0..1,true] = "; p na[0..1,true]
  15|print "na[true,-1]";      p na[true,-1]
  16|print "na[[0,2,0],true]"; p na[[0,2,0],true]
  17|
  18|include NMath
  19|print "sqrt(na) = "; p sqrt(na)
  20|print "log10(na + 1) = "; p log10(na + 1)

実行結果は下のようになる。

 % ruby na1.rb
 na = NArray.float(3,2):
 [ [ 0.0, 1.0, 2.0 ],
   [ 3.0, 4.0, 5.0 ] ]
 nb = NArray.float(1,2):
 [ [ 1.0 ],
   [ 10.0 ] ]
 nc = NArray.sfloat(12):
 [ 0.09706, 0.03895, 0.3794, 0.6282, 0.6721, 0.3771, 0.2122, 0.8495, ... ]
 na + nb = NArray.float(3,2):
 [ [ 1.0, 2.0, 3.0 ],
   [ 13.0, 14.0, 15.0 ] ]
 na >= nb = NArray.byte(3,2):
 [ [ 0, 1, 1 ],
   [ 0, 0, 0 ] ]
 na[na>=nb]=999 --> na = NArray.float(3,2):
 [ [ 0.0, 999.0, 999.0 ],
   [ 3.0, 4.0, 5.0 ] ]
 na[0..1,true] = NArray.float(2,2):
 [ [ 0.0, 999.0 ],
   [ 3.0, 4.0 ] ]
 na[true,-1]NArray.float(3):
 [ 3.0, 4.0, 5.0 ]
 na[[0,2,0],true]NArray.float(3,2):
 [ [ 0.0, 999.0, 0.0 ],
   [ 3.0, 5.0, 3.0 ] ]
 sqrt(na) = NArray.float(3,2):
 [ [ 0.0, 31.607, 31.607 ],
   [ 1.73205, 2.0, 2.23607 ] ]
 log10(na + 1) = NArray.float(3,2):
 [ [ 0.0, 3.0, 3.0 ],
   [ 0.60206, 0.69897, 0.778151 ] ]

プログラムの第 2-4 行目では、3つの NArray を初期化する。3 行目の NArray.to_na メソッドは、このように配列を NArray に変換するだけなく、 多彩な機能を持つ。例えば、バイナリーデータを入れた String と、配列の取 るべき形を与えて初期化することも出来るため、次のように、バイナリーファ イルからデータを読むこともできる:

 file = File.open('hoge.bin', 'rb')
 s = file.read( 4*nx*ny )
 na = NArray.to_na(s, NArray::SFLOAT, nx, ny)

バイトオーダーの変換も出来る:

 na = na.ntoh    # ネットワークバイトオーダー (big endian) から、その
                 # 計算機のエンディアンに変換。その他変換関数多数あり。

さて、上のプログラムの 9 行目では足し算、10 行目では大小比較を行っている。 2項演算は、同じ形の配列に対する要素毎の演算となる。しかし、長さ1の次 元は任意の長さの次元に拡張されることになっているため、上のように 3×2 の配列と 1×2 の配列の演算が行えるのである。なお、真 / 偽は 1 / 0 の整 数で表される (バイト型以外も可)。これをマスクとして、部分配列に代入する こともできる (11 行目)。もちろん、部分配列の指定は、14-16 行目のように添 字でも行える。from..to のように Range オブジェクトで範囲を指定でき、 true は全選択、すなわち 0..-1 と等価である。

三角関数等は、NArray 対応した Math モジュールである NMath に収められて いる。18-20 行目はその使用例である。NArray には、他にも数値型多次元配 列にふさわしいメソッドが多数用意されている。例をいくつか挙げる:

 na.mean(dim,..)    # 指定した次元に関する平均. 引数なしなら全要素の平均.
 na.stddev(dim,..)  # 同標準偏差 (分母に n-1 を使う不偏分散ベース)
 na.transpose(dim0, dim1,..)  # 次元の入れ替え
 na.reshape!(size0,..)  # 配列の形の変更 (全要素数は不変).

などなど。ここで次元は、最初 (左端) の次元を 0 とする整数で指定する。

なお、NArray では、部分配列の指定においてステップを指定して間引くこと が出来ないというのが、多次元配列クラスとして唯一見劣りすると筆者が思う 点である。実は拡張ライブラリー内には間引き機能も実装されているのだが、 それを Ruby で表現する美しい方法がなかなか見当たらないというのが理由だ そうである (例えばステップ付きの Range クラスのようなものがあれば解決す るのだが)。今後に期待したい。現状では、飛び飛びの添字の列を与えるとい うことで対応することになる:

na2.rb

   1|require "narray"
   2|na = NArray.int(10).indgen!
   3|str = 1
   4|stop = 7
   5|step = 2
   6|p na[ str + step * NArray.to_na([0..((stop-str)/step)]) ]

指定した次元に関する操作

上記の mean など、特定の次元に対する操作は、数値計算ではよく使われる。 Ruby では [ ] がメソッドであり、さらに配列を '*' で展開して引数にでき るので、汎用性のあるメソッドが簡単に実現できる。次のプログラムは指定し た次元に沿っての差分を計算する。

na3.rb

   1|require "narray"
   2|
   3|class NArray
   4|  # differentiate along the dim-th dimension:
   5|  #   self[.., 1..-1, ..] - self[.., 0..-2, ..]
   6|  def diff(dim)
   7|    idx1 = [true]*dim + [1..-1, false]
   8|    idx2 = [true]*dim + [0..-2, false]
   9|    self[*idx1].sbt!( self[*idx2] )   # sbt! : destructive subtraction
  10|  end
  11|end
  12|
  13|nx = 5
  14|ny = 4
  15|x = NArray.float(nx,1).indgen!
  16|y = NArray.float(1,ny).indgen!
  17|na = x * y
  18|p na, na.diff(1)

7、8 行目の false は、0 以上の任意個の true を表す「ゴム次元」である。 実際何個に対応するかは、行列の次元数に依存する。つまり、7 行目は

 idx1 = [true]*dim + [1..-1] + [true]*(self.rank-dim-1)

と同等である。なお、ゴム次元をもつ配列を、文法としてサポートした言語に は yorick というものがある(ゴム(rubber)次元というのも、yorick の用語 を借用して筆者が勝手にしている呼び方である)。

NArray は [ ] で部分配列を切り出すと独立した新しいオブジェクトを作る。 このため、9 行目の括弧内の self[*idx2] で新規にメモリー割り当てが行われ、 その後はごみとして GC に委ねられるという無駄が生ずる。もしこれが良く使 うメソッドで効率を求めるという話であれば、C で書き直すというのが一法で ある。いずれにしろ、数値微分に使うためには、境界条件の処理も考える必要 があるが。

上のプログラムの実行結果を示す。次元1、すなわち第2次元目に関する差分が 計算できていることが分かるであろう。

 % ruby na1.rb
 NArray.float(5,4):
 [ [ 0.0, 0.0, 0.0, 0.0, 0.0 ],
   [ 0.0, 1.0, 2.0, 3.0, 4.0 ],
   [ 0.0, 2.0, 4.0, 6.0, 8.0 ],
   [ 0.0, 3.0, 6.0, 9.0, 12.0 ] ]
 NArray.float(5,3):
 [ [ 0.0, 1.0, 2.0, 3.0, 4.0 ],
   [ 0.0, 1.0, 2.0, 3.0, 4.0 ],
   [ 0.0, 1.0, 2.0, 3.0, 4.0 ] ]

行列計算

最後に、行列計算の例を示す。NArray のサブクラス NMatrix には、ピボット 選択付の LU 分解が実装されており、たちの良い線型方程式の解や逆行列を計 算することができる。下で題材に使った Hilbert 行列というのは、サイズが 大きくなるにつれ急速に特異行列に近づくため、非常に解きにくいことで有名 である。以下では、誤差の目安として、逆行列と元の行列の積と、単位行列の 差の絶対値の最大値のみを示することにする。

nmatrix_lu.rb

   1|require "narray"
   2|
   3|def hilbert(n)
   4|  h = NMatrix.float(n,n)
   5|  for i in 0...n
   6|    for j in 0...n
   7|      h[j,i] = 1.0 / (i+j+1)
   8|    end
   9|  end
  10|  h
  11|end
  12|
  13|[3,6,12,13].each do |n|
  14|  begin
  15|    h = hilbert(n)
  16|    e = NMatrix.float(n,n).unit
  17|    hi = e / h
  18|    print "n=#{n}. Error in h*h^-1 : ", (e - h*hi).abs.max, "\n"
  19|  rescue
  20|    print "n=#{n}. COMPUTATION FAILED: ",$!,"\n"
  21|  end
  22|end

実行結果は下記。倍精度計算であるが、n=12 では誤差がそれなりに大きく n=13 では特異ベクトルと判別されて計算できてない。これは、そのような数 値的に難しい問題だからである。

 % ruby nmatrix_lu.rb
 n=3. Error in h*h^-1 : 3.68455266297474e-15
 n=6. Error in h*h^-1 : 1.40474298859772e-10
 n=12. Error in h*h^-1 : 0.0399007797241211
 n=13. COMPUTATION FAILED: singular matrix, status=1

なお、NArray は Ruby オブジェクトを要素とすることが出来るので、 今回の問題では有理数を使えば正確に計算することができる。 ソースはこちら:nmatrix_lu_rational.rb。 実行結果は、下記のように誤差ゼロである。全ての計算が Ruby のメソッドを呼び出す形で行われるので時間はかかる。 と言っても、筆者のパソコンでの実行時間は約1秒である。

 % ruby nmatrix_lu_rational.rb
 n=3. Error in h*h^-1 : 0
      max and min of h^-1 : 192, -180
 n=6. Error in h*h^-1 : 0
      max and min of h^-1 : 4410000, -3969000
 n=12. Error in h*h^-1 : 0
      max and min of h^-1 : 3659449159080000, -3110531785218000
 n=13. Error in h*h^-1 : 0
      max and min of h^-1 : 106518477825760000, -103492384705710000

逆行列の要素の値は非常に大きく、n=12 でも Bignum (任意長整数) の領域に 入っている。言い換えれば、Bignum をサポートする Ruby ならではの妙技で ある。

レポート RubyDCL

登録データ

概要

DCL (Dennou Club Library) は、FORTRAN77 で書かれた可視化ライブラリーである。RubyDCL は、移植性など を考慮して C 化した、DCL の C バージョンを使うラッパーであり、 オリジナルと一対一に対応する。引数の仕様も基本的に同一で あるが、配列サイズは引数にせず、出力はすべて返り値で受け、作業領域は内 部で生成するという点が異なる。マニュアルも、FORTRAN 版から機械的に生成 されたものである。

RubyDCL は、配列データとして NArray、Array 及び NArrayMiss を受け入 れる。NArrayMiss は、NArray を用いて作成された、データ欠損のあり得る配 列クラスである。平均は欠損をよけて行うなどの特徴がある。なお、DCL 内で 欠損値処理させるだけであれば、NArray のみで十分である。

DCL がサポートする可視化法は、折れ線やヒストグラム、等高線、色塗り、ベ クトル図などの1、2次元のグラフィックである。それらを 3次元空間に投影す ることもできるが、等値面作成やボリュームレンダリングといった本格的な3 次元可視化機能はない。その分、1、2次元可視化は多彩で、出版に耐えるきめ 細かなレイアウトが可能である。なお、一般的な数値データの可視化機能に加 え、データの欠損に対応でき、さらに地図投影など地球科学固有の問題に対処 する。

RubyDCL を実際に使いこなそうと思うと、慣れるまではマニュアルと格闘する 必要があるであろう。それが辛いところではある。ただ、ホームページ上に、 充実したチュートリアル「ごくらくDCL」、「らくらくDCL」が置いてあり、配 布パッケージには、両チュートリアルの全プログラムを含む豊富なデモプログ ラムが含まれるので、参考になるであろう。今後、チュートリアルの実行結果 の図から、描画プログラムを「逆引き」する web ページも用意する計画があ る。

サンプル

1次元描画

まずは縦軸、横軸の値を 1次元配列で与える折れ線プロットの例を示す。

dcl_hop.rb

   1|require "narray"
   2|require "numru/dcl"
   3|include NumRu                 # Most dennou-ruby products are in it
   4|include NMath                 # Math module for NArray
   5|
   6|n = 400
   7|dt = 2*PI/(n-1)
   8|t = NArray.sfloat(n).indgen! * dt        # 0 <= t <= 2*pi
   9|x = 1e2 *sin(4*t)
  10|y = 1e-3*cos(5*t)+6
  11|
  12|DCL::gropn(1)                      # 1=>Display; 2=>PS file; 4=>Gtk
  13|DCL::grfrm                               # initilize a frame (page)
  14|DCL::ussttl('X-TITLE', 'x-unit', 'Y-TITLE', 'y-unit')
  15|DCL::usgrph(x, y)
  16|DCL::grcls

これを実行すると下に示す図が表示される。

電脳 Ruby プロジェクト製品の多くは、名前空間確保のため NumRu (語源は Numerical Ruby) というモジュールに入れてある。利用の際には 2 行目のよ うに require "numru/dcl" し、特に不都合がなければ 3 行目のように include NumRu して、むき出しにすれば良かろう。DCL のメソッドはすべてNumRu::DCL というモジュールに入っている。4 行目の NMath は NArray 対応の Math モ ジュールである(前述)。

6-10 行目ではプロットすべきデータを用意してる。t を媒介変数として x, y がいわゆるリサジュー図形を表す。y は 6 という値を中心に振幅 0.001 で振 動するという、座標軸を描く上ではいやらしいデータである。しかし、DCL は 図に示すように、うまく対処する(図の左上の括弧内を見よ)。

12-16 行目が本題の可視化部分である。最初に (12 行目) グラフィックデバイス を初期化する。引数 1 は画面上での表示を指示する。多くの OS では X window システムを使うが、Windows 版では Windows ネーティブな描画になる。 引数が 4 の場合も同様に画面に出すが、Gtk を用いる (Gtkのないシステムで は 4番は利用できない)。1, 4 ともに画像ダンプが可能である。下のスクリー ンショットと違い、ウィンドー枠は含まない. 引数が 2 の場合 (Windows 版 では 3 の場合) は Postscript ファイルを出力する。13 行目で作図するペー ジを初期化し、15 行目で実際の作図が行われる。変数 x, y をそれぞれ横軸、 縦軸にとる折れ線を描き、座標軸を描いた上で、折れ線で描画する。その際、 14 行目のように、次の作図で使う縦・横軸のタイトル及び単位を設定するこ とができる。16 行目はグラフィックデバイスの終了処理である。

dcl_hop_sml.png
dcl_hop.rb 実行結果のスクリーンショット

以上、初期化と終了処理が必要なのがやや煩わしいかもしれないが、とりあえ ずお任せで図を描くだけであれば、まずまず簡単にできることがわかる。さら に、線種を変えたり、描画範囲を陽に設定したり、座標軸をカスタマイズした り、1ページに複数の図を出したりといろいろ出来るが、マニュアルを見るか デモプログラムを参照しなければならないのは、プログラミングで図を描く場 合の辛いところである。反面、一度描いた図と似たような図の作成は、プログ ラムの改変で行えるので、比較的楽である。

2次元描画

今度は、2次元配列を使って、等値線 (コンター) 描きと色の塗り分けをしてみ よう。描画範囲も陽に設定するなど、上の例よりちょっと複雑なプログラムに なる。

dcl_contour.rb

   1|=begin
   2|=dcl_contour.rb: 等高線と色塗りのサンプル
   3|使用法
   4| % ruby dcl_contour.rb [-ps] [-dump]
   5|履歴
   6| * 水田亮 作成 (「ごくらくDCL」デモプログラムを改変)
   7| * 堀之内 改造 2005/04
   8|=end
   9|
  10|require "numru/dcl"
  11|include NumRu
  12|include NMath
  13|
  14|nt = 50
  15|nz = 50
  16|tmin, tmax = 0.0, 5.0
  17|zmin, zmax = 20.0,50.0
  18|
  19|t = NArray.sfloat(nt+1,    1).indgen! * (tmax-tmin)/nt    # 横軸
  20|z = NArray.sfloat(   1, nz+1).indgen! * (zmax-zmin)/nz    # 縦軸
  21|
  22|uz = exp(-0.2*z)*sqrt(z)
  23|tz = -2.0*exp(-0.1*z)
  24|u = uz*sin(3.0*(tz+t))               # (nt+1,nz+1)の配列になる。
  25|
  26|if ARGV.index("-ps")
  27|  # コマンドラインオプションに-psが含まれていればPSファイルを作る
  28|  DCL::gropn(2)     # Windows では 3 にすること
  29|else
  30|  # ウィンドウを開く
  31|  DCL::swpset('iheight',700)           # 画面の縦方向のピクセル数指定
  32|  DCL::swpset('iwidth',700)            # 画面の横方向のピクセル数指定
  33|  DCL::gropn(1)     
  34|  DCL::swpset('ldump', true) if ARGV.index("-dump")   # 画面のダンプ
  35|end
  36|
  37|DCL::grfrm
  38|DCL::grswnd(tmin, tmax, zmin, zmax)         # 左・右・下・上の座標値
  39|DCL::uspfit                                 # あとはおまかせ
  40|DCL::grstrf                                 # 座標の確定
  41|
  42|DCL::ussttl('TIME', 'YEAR', 'HEIGHT', 'km') # 軸のタイトルと単位の設定
  43|DCL::uelset("ltone",true)                   # 色つきにする
  44|DCL::uetone(u)                              # ぬりわけ
  45|DCL::usdaxs                                 # 座標軸を描く
  46|DCL::udcntr(u)                              # 等値線をひく
  47|DCL::uxmttl('t','Idealized QBO',0.0)        # title. top に中央合せ(0.0)
  48|
  49|DCL::grcls

19-24 行目では、表示すべきデータを生成している。やや長いが、ソースがほ ぼ数式どおりに書かれているので分かりやすい。NArray のレポートで説明し たように、長さ 1 の次元は2項演算時に任意の長さに拡張されることを利用し ていることに注意。

26-35 行目ではグラフィックデバイスの初期化を行っている。オプションに応 じて PS ファイルに出力したり、画像ダンプを取るようにしている。37-40 行 目では、ページを初期化し (37 行目)、画面の4隅の座標値を設定し (38 行目)、 それを画面上のどの範囲に表示するかはお任せとし (39 行目)、座標系の設定を 確定させる (40 行目)。42-47 行目が実際の描画である。まず軸のタイトルを設 定し (42 行目)、次に行う塗り分けをカラーですることを指定し (43 行目)、塗り 分けし (44 行目)、座標軸を描画し (45 行目)、等高線を書き (46 行目)、図の上に 大きめのタイトルを書く (47 行目)。以上で現れる、ページの初期化や座標軸の 描画命令の多くが、上の折れ線図の例には現れなかったのは、usgrph という 複合メソッド内で自動的に呼ばれていたからである。

dcl_contour_sml.png
dcl_contour.rb 実行結果のスクリーンショット

データ欠損の処理、地図投影、3次元、etc etc

RubyDCL 及び DCL のチュートリアルページ「ごくらくDCL」、「らくらくDCL」 を参照のこと。例えばこんな図が描ける:

dcl_map_proj.png dcl_3D.png
この図を作ったプログラム: dcl_map_proj.rb この図を作ったプログラム: dcl_3D.rb

その他の可視化ライブラリー

Ruby-VTK

C++ ベースのフリーな 3次元可視化ライブラリー Visulaization Tool Kit (VTK; ホームページは http://public.kitware.com/VTK/) の Ruby インター フェースである。 SWIG による自動生成で、VTK の全機能が利用可能 である。本家の VTK は、TK を用いてグラフィカルユーザーインターフェース

(GUI) が作れるが、Ruby VTK では Ruby/Gtk で GUI を作れる (Ruby/Tk も使

えるはずである)。

下の図は、Ruby-VTK 作者が書いたサンプルプログラムを使って、ある数値デー タの等値面を描画した際のスクリーンショットである。 vtk_surface_sml.png

Ruby/PGPLOT

グラフィックライブラリー PGPLOT のラッパー。NArray の作者により作られ たもので、データは NArray で与える。機能的には DCL と同様なものである ようだが、筆者は使ったことがないので、間違った紹介をしないよう、名前を 挙げるに留める。

その他

ヒストグラムや円グラフなどが美しく書ける gdchart のラッパーなど、複数 のライブラリーが RAA の Library/Graphics カテゴリーに複数登録されて いる。他に、パイプラインを使って gnuplot を呼ぶ Ruby Gnuplot もある。 次に紹介する数値計算ライブラリー Ruby/GSL には、GNU plotutils を使った 可視化機能がある。

レポート Ruby/GSL

登録データ

概要

GNU Scientific Library (GSL) は、C で書かれた、フリーな数値計算ライブ ラリーである。様々な数値計算が実装され、現在も開発が行われている。 Ruby/GSL はそのほとんどカバーする Ruby インターフェースである。GSL の 諸関数の純粋なラッパーに留まらず、使い勝手と CPU 効率を考えた拡張がな されている優れたライブラリーである。例えば、本家の GSL では、補間を行 う関数は 1点のみに対して行う。Ruby/GSL でも同じことが出来るが、さらに、 NArray 等の配列を与えると全要素を順に処理し、結果をまとめて同じ型の配 列で返すというようになっている。

GSL 本家は、データ保持に構造体を多用する。例えば、行列は gsl_matrix、 ベクトルは gsl_vector, 複素数は gsl_complex, 特殊関数は gsl_sf_result と言った具合いで、例えば行列の構造体ならサイズも込になっているなど、 stand alone なデータ構造になっている。Ruby/GSL では、これらを GSL::Matrix, GSL::Vector, GSL::Complex, GSL::Sf::Result などとクラス化 する。一方で、GSL では、1次元配列は剥き出しのポインターと長さの組で扱 われる場合も多い。Rub/GSL は、このような場合は GSL::Vector で対応する。

Ruby/GSL は、GSL の構造体をラップしたクラスだけでなく、NArray も利用し やすいように出来ている。NArray と GSL::Matrix 等の間の相互変換をサポー トするのみならず、多くの関数が、引数に直接 NArray を取ることが出来るよ うになっている。その際、可能な限り、データポインターをコピーなしに構造 体間で受け渡しして、計算効率が落ちないように配慮されている。

GSL 本体にはない Ruby/GSL の独自機能として、GNU graph (GNU plotutils に含まれる描画コマンド) を用いるグラフィックスがある。折れ線などによる プロットやヒストグラム、関数の表示などができる。お任せのレイアウトであ れば、GSL::Vector などに用意された graph メソッドを呼ぶだけで良く、さ らに graph のオプションを与えてカスタマイズできる。Ruby/GSL のホームページにサンプルとスクリーンショットがあるので参照されたい。

サンプル

GSL の機能を挙げ出したら切りがない。詳しくはホームページを見て頂くこと にして、ここでは一例のみ、NArray の紹介で取り上げた LU 分解による逆行 列を、GSL を使って求めてみることにする。GSL の構造体を使う場合と、 NArray を使用する場合の二通りを示す。

なお、GSL のドキュメントでは、GSL の行列計算は規模の小さい場合向きであ り、大規模計算には LAPACK を使うと良いとされている。このような場合は、 後述の Linalg (LAPACKベース) や Ruby-SSL2 を使うと良いだろう。

GSL のみ使う場合

次のプログラムは、前述の nmatrix_lu.rb を NArray を使わずに書き直した ものである。

gsl_lu.rb

   1|require "gsl"
   2|
   3|def hilbert(n)
   4|  # 実は GSL::Matrix にはクラスメソッドとして hilbert があるのだがデモのため
   5|  h = GSL::Matrix.new(n,n)
   6|  for i in 0...n
   7|    for j in 0...n
   8|      h[i,j] = 1.0 / (i+j+1)
   9|    end
  10|  end
  11|  h
  12|end
  13|
  14|[3,6,12,13].each do |n|
  15|  h = hilbert(n)
  16|  hi = GSL::Linalg::LU.invert(h)
  17|  e = GSL::Matrix.unit(n)
  18|  diff = e - h*hi
  19|  err = [ diff.max, -diff.min ].max      # diff.abs.max if rbgsl >= 1.6.3
  20|  print "n=#{n}. Error in h*h^-1 : ", err, "\n"
  21|end

GSL::Matrix の値の代入は NArray と同じようにできることが分かる。ただし、 こちらは、添字は普通に行、列の順である。例で用いたのは対称行列なので、 間違えても結果は変わらないが。Ruby/GSL のメソッドとしては、LU 分解だけ 行うものから、線形方程式を解くもの、そしてここで用いた逆行列を求めるメ ソッドと、各種取りそろえているのは、NArray の LU 分解と同様である。 ただ、NArray と違って、GSL では線形方程式の解法だけでも多種あり、他に も固有値計算や特異値分解などがある。

Ruby/GSL 1.6.2 までは GSL::Matrix の各要素の絶対値を取るメソッドがない ので、誤差評価は 19, 20 行目のようにした。実行結果は 下記。GSL ではより大きなサイズ (n=13) でも解は求まってはいるが、精度が 悪いので、今のケースでは NArray の LU 分解とほぼ同等な結果が得られたと 思って良い。

 % ruby gsl_lu.rb
 n=3. Error in h*h^-1 : 3.33066907387547e-15
 n=6. Error in h*h^-1 : 5.9657168094418e-11
 n=12. Error in h*h^-1 : 0.0218987464904785
 n=13. Error in h*h^-1 : 0.346458435058594

GSL と NArray の連携

Ruby/GSL の前に NArray をインストールしておくと、Ruby/GSL でも NArray を使うことができる。次に出る 1.6.3 以降では、自動的に NArray を検出す るが、検出されない場合や 1.6.2 まででは、

 % ruby setup.rb config -- --with-narray-include=<directory where 'narray.h' exists>

という形で configure する。

NArray 対応の Ruby/GSL では、行列に関し以下のメソッドが定義されている。

  GSL::Matrix#to_na          NArray へ変換
  GSL::Matrix#to_na_ref      参照タイプの NArray へ変換 (データを共有)
  NArray#to_gm               GSL::Matrix へ変換
  NArray#to_gm_view          参照タイプの GSL::Matrix へ変換 (データを共有)

上で ref と view という異なる用語が使われているのは、NArray と GSL そ れぞれの用法による。さらに、1.6.3 以降では以下も使えるようになる。

  GSL::Matrix#to_nm          NMatrix へ変換
  GSL::Matrix#to_nm_ref      参照タイプの NMatrix へ変換 (データを共有)

なお、NMatrix は NArray のサブクラスなので、現状でも to_gm は使える。

それでは、上の gsl_lu.rb を NArray (NMatrix) を使う形で書き直したもの を示す。

gsl_lu_narray.rb

   1|require "narray"
   2|require "gsl"
   3|
   4|def hilbert(n)
   5|  h = NMatrix.float(n,n)
   6|  for i in 0...n
   7|    for j in 0...n
   8|      h[j,i] = 1.0 / (i+j+1)
   9|    end
  10|  end
  11|  h
  12|end
  13|
  14|[3,6,12,13].each do |n|
  15|  h = hilbert(n)
  16|  ghi = GSL::Linalg::LU.invert(h.to_gm_view)
  17|  hi = NMatrix.refer( ghi.to_na )       # Simply ghi.to_nm if RubyGSL >= 1.6.3
  18|  e = NMatrix.float(n,n).unit
  19|  print "n=#{n}. Error in h*h^-1 : ", (e - h*hi).abs.max, "\n"
  20|end

ここでは NArray と GSL::Matrix の間の変換は陽に行った。しかし、NArray を直接与えて動作するメソッドも多く、かつ増え続けている。その場合、 入力引数が NArray なら、出力も NArray になるよう工夫されている。

その他の数値計算ライブラリー

Ruby-SSL2

Fortran ベースの商用数値計算ライブラリー、富士通 SSL2 の Ruby インター フェースであり、SSL2 の全機能が利用できる。SSL2 は、スーパーコンピュター で鍛えられてきた歴史ある優れた数値計算ライブラリーで、現状では性能は一 般に GSL より勝っているであろう。SSL2 は富士通 の Fortran コンパイラー に付属し、マニュアルも充実している。なお、Ruby-SSL2 では、配列には NArray を利用する。組み込みの Array も利用可能である。

Linalg

有名で実績のある線形代数ライブラリー LAPACK を用いるためのライブラリー。LAPACK のラッパーと、それをもとにした行列計算ライ ブラリーからなる。

Linalg の LAPACK モジュールは、FORTRAN77 による LAPACK のオリジナルライ ブラリーを1対1にラップする。その戦略は、Ruby/GSL や Ruby-SSL2 のそれ とは違い、ポインターのアドレスを Ruby の Integer オブジェクトとしてラッ プするというものである。FORTRAN77 の変数は C から見れば全てポインターで ある。このため、各引数のアドレスを Integer オブジェクトにすることで、 Ruby と FORTRAN 副プログラムの橋渡しをすることができるのである。これを 利用して、Linalg では簡単なフィルターで LAPACK のヘッダーファイルを処 理するだけで、LAPACK の全機能搭載を実現している (LAPACK は C 版もある のでヘッダーファイルもあるのである)。

この LAPACK モジュールは、そのままでは使いにくいし、危険でもあるので、 その上にラッパーを被せてある. それが Linalg の上位ライブラリーである。 例えば行列に関しては、[行サイズ, 列サイズ, 中身のポインター] からなる "?Matrix" (?=D,C,etc) クラスが用意されている. それに定義されているデー タ領域のポインターを返すメソッドを利用して、Ruby レベルで LAPACK モジュー ルのメソッドを呼ぶ。なお、上位ライブラリーにおいては、今のところ LAPACK の全サブルーチンまではサポートされてない。

上述のように LAPACK モジュールの関数はポインターベースであるため、 NArray / NMatrix や GSL::Matrix などといった、他の拡張ライブラリーで使 うことも可能である。ただし、LAPACK のデータ構造は column-major である ため、行列についてはデータの詰め替え(transpose)が必要であり、また、デー タ領域のポインターを Integer で返すメソッドを追加する必要がある. NArray で LAPACK がフルに使える、使いやすいモジュールが欲しいところで ある。

レポート RubyNetCDF

登録データ

概要

NetCDF は、気象・海洋分野でよく用いられるバイナリー数値データのフォー マットであり、そのための入出力ライブラリーである。NetCDF は名前ベース で内容にタグづけしており、実際のバイナリー構造は隠蔽される。よって、専 用の入出力ライブラリーを使うことが前提となっている。ライブラリーは C で書かれており、その上にいくつもの言語へのバインディングが行われている。 一方、独立に Java による入出力ライブラリーも提供されている。

RubyNetCDF は、C 版の NetCDF ライブラリーを用いた拡張ライブラリーであ る。単なるラッパーではなく、NetCDF のデータ構造を考えてクラスライブラ リーを構成している。NetCDF はオブジェクト指向に適したデータモデルを取っ ているので、Ruby 版の構成はごく自然な形で行われている。次にそれを説明 する。

NetCDF のデータモデル

NetCDF の構成を示す。[ ] 内は RubyNetCDF におけるクラス名である。 NetCDFClassDiagram.png データセット中には、0 次元 (即ちスカラー) またはそれ以上の次元の数値配列 である「変数」が定義できる。変数同志は名前で区別される。データセット及 び変数は、任意個の「属性」を持てる。属性とは名前と値の組であり。値は文 字列、または同一型で任意個の数値である。変数の各次元の長さは、数値で指 定するのでなく、名前付の「次元」で定義される。「次元」はそれぞれ長さを 持つが、可変な「無制限次元」も定義できる。

以上は、バイナリー数値データを持つための一般的な枠組みである。NetCDF の特徴は、この上に「物理」的なデータを保持するための仕組みを定めている ことにある。それを使わなくてもいいという柔軟性がある一方、作法を強制で きないのは弱点とも言える。ともあれ、仕組みはこうである:

  • 座標の関連づけ
    • 次元名と同じ名前の一次元変数があれば、座標値を表すものとする。例えば、経度に関し 0,10,20,..350 とサンプリングしたデータに対し、"lon" という名の長さ 36 の「次元」を定義し、同名の1次元配列に[0,10,..] の値を収める。物理量を収める配列の定義にも次元 "lon" を用いることで、物理量値の変数と座標変数が関連づけられる。''
  • 属性の標準化
    • NetCDF の属性は任意に付けることが出来るが、単位を "units" で表すなどといった標準化が行われており、実際に多くの NetCDF 作者が従っている。このため、機械的な自動処理で、基本的な情報が引出せる。これを基本に、さらに複雑なデータに対応するための標準化も行われている。

サンプル

データ書出し

NetCDF ファイルをゼロから作ってみよう。その前に、NetCDF ファイルの中身 を調べるコマンドを紹介する。下記の ncdump というのは、(本家の) NetCDF と共にインストールされるコマンドであり、バイナリーの NetCDF を、 CDL (Common Data Language) というテキスト表現に変換する。順序は逆にな るが、これから作るファイルを ncdump したものを示す。但し、ファイルの構 成のみを示すため -h オプションをつける (ヘッダー表示モード):

 % ncdump -h test.nc
 netcdf test {
 dimensions:
         x = 50 ;
         y = 25 ;

 variables:
         float x(x) ;
                 x:long_name = "longitude" ;
                 x:units = "degrees" ;
         float y(y) ;
                 y:long_name = "latitude" ;
                 y:units = "degrees" ;
         float u(y, x) ;
                 u:long_name = "wind speed" ;
                 u:units = "m/s" ;

 // global attributes:
                 :history = "Created by netcdf_create.rb  2005-04-25" ;
 }

これを作ったプログラムを示す。

netcdf_create.rb

   1|require "numru/netcdf"
   2|include NumRu
   3|include NMath
   4|
   5|#   -- create a new dataset --
   6|file = NetCDF.create("test.nc")             # .nc is the standard suffix
   7|
   8|#   -- define data --
   9|nx, ny = 50, 25
  10|xdim = file.def_dim("x",nx)
  11|ydim = file.def_dim("y",ny)
  12|require "date"
  13|file.put_att("history","Created by #{$0}  #{Date.today}")
  14|
  15|xvar = file.def_var("x","sfloat",[xdim])        # x coordinate variable
  16|xvar.put_att("long_name","longitude")
  17|xvar.put_att("units","degrees")
  18|
  19|yvar = file.def_var("y","sfloat",[ydim])        # y coordinate variable
  20|yvar.put_att("long_name","latitude")
  21|yvar.put_att("units","degrees")
  22|
  23|uvar = file.def_var("u","sfloat",[xdim,ydim])   # 2D along x and y 
  24|uvar.put_att("long_name","wind speed")
  25|uvar.put_att("units","m/s")
  26|
  27|file.enddef                                  # end the define mode
  28|
  29|#   -- create data values for test --
  30|x = NArray.sfloat(nx, 1).indgen! * (360.0/(nx-1))
  31|y = NArray.sfloat(1, ny).indgen! * (180.0/(ny-1)) - 90.0
  32|x_rad = x * (PI/180)
  33|y_rad = y * (PI/180)
  34|u = (30 + 5*cos(2*x_rad)) * ( 0.2 - cos(3*y_rad) )
  35|
  36|#   -- write data values --
  37|xvar.put( x )
  38|yvar.put( y )
  39|uvar.put( u )
  40|
  41|file.close

やや長いが、読めば何をしてるかは分かるであろう。 なお、2 行目の NumRu については、RubyDCL の節を、 3 行目の NMath については NArray の節を参照のこと。

イテレーター

RubyNetCDF の各クラスには Ruby ならではのイテレーターが定義されている。 下のプログラムは、任意の NetCDF ファイルのコピーを作る。 NetCDF のデータモデル と比較すると容易に理解できるだろう。

netcdf_copy.rb

   1|#  Makes a copy of a NetCDF file:
   2|#  % ruby netcdf_copy.rb path_from path_to
   3|
   4|require "numru/netcdf"
   5|include NumRu
   6|usage = "\n\nUSAGE:\n% ruby #{$0} path_from path_to\n"
   7|raise usage if ARGV.length != 2
   8|path_from, path_to = ARGV
   9|from = NetCDF.open(path_from)
  10|to = NetCDF.create(path_to)
  11|from.each_dim{|dim| to.def_dim(dim.name, dim.length_ul0)} # len==0 if unlimited
  12|from.each_att{|att| to.put_att(att.name, att.get)}
  13|from.each_var do |var|
  14|  newvar = to.def_var(var.name, var.ntype, var.dim_names)
  15|  var.each_att{|att| newvar.put_att(att.name, att.get)}
  16|end
  17|to.enddef
  18|from.each_var{|var| to.var(var.name).put(var.get)}
  19|to.close

ちなみに、NetCDF の C ライブラリーは、変数の型毎に関数が用意されている ため、上のプログラムを C で書くと条件分岐の嵐となり、非常に長いプログ ラムとなる。

読み込み

次は、上の netcdf_create.rb により作成したファイルを読み込み、図示する。 可視化には RubyDCL を用いる。

netcdf_plot.rb

   1|require "numru/netcdf"
   2|require "numru/dcl"
   3|include NumRu
   4|
   5|#  -- read data --
   6|file = NetCDF.open("test.nc")          # made with netcdf_create.rb
   7|xvar = file.var("x")                   # get a NetCDFVar object named "x"
   8|yvar = file.var("y")
   9|uvar = file.var("u")
  10|x = xvar.get                           # read values of the variable
  11|y = yvar.get
  12|u = uvar.get
  13|xname = xvar.att('long_name').get      # get a NetCDFAtt and read values
  14|yname = yvar.att('long_name').get
  15|uname = uvar.att('long_name').get
  16|xunits = xvar.att('units').get
  17|yunits = yvar.att('units').get
  18|uunits = uvar.att('units').get
  19|
  20|# -- graphics --
  21|DCL.gropn(1)
  22|DCL.sgpset('lfull',true)       # use the whole window even if not a square
  23|DCL.uzfact(0.8)                             # --> a bit smaller characters
  24|DCL.grfrm
  25|DCL.grswnd(x.min, x.max, y.min, y.max)        # physical boundaries
  26|DCL.grsvpt(0.2, 0.8, 0.2, 0.6)                # viewport on window
  27|DCL.grstrf                                    # fix the "transform"
  28|DCL.uelset("ltone",true)                      # color shading from now on
  29|DCL.uetone(u)                                 # shading
  30|DCL.ussttl(xname, xunits, yname, yunits)      # set axis titles
  31|DCL.usdaxs                                    # axis
  32|DCL.udcntz(u)                                 # contour
  33|DCL.uxmttl('t', "#{uname} [#{uunits}]", 0.0)  # title at the top
  34|DCL.grcls

これを実行すると、次のような図が表示される。

netcdf_plot_sml.png
netcdf_plot.rb 実行結果のスクリーンショット

さて、このプログラムでは変数 u, x, y を陽に指定して読み込んでいる。し かし、次元名=変数名な変数を座標変数とする規約から、u さえ指定すれば、 x, y は自動的に手繰れるはずである。なのに別々に扱うのは冗長ではないか?  そこを自動化するライブラリーを作れば、図を書くのも、そしてデータの解析 ももっと楽になるはずである。次に紹介する GPhys はそれを実現したもので ある。

レポート GPhys

登録データ

概要

GPhys は、格子点上の物理 量のためのクラスライブラリーである。入出力から、数値解析、可視化まで幅広 く扱う包括的ライブラリーである。 但し、作者の都合で、実際には大気・海洋の研究に使われている。利用者は限ら れるであろうが、これまで紹介してきたライブラリーをもとに、便利なライブラ リーを開発する例として、簡単に紹介する。

GPhys で基本になるのは、RubyNetCDF の項で示したようなデータモデルであ る。すなわち、各物理量は、データの値だけでなく、座標や単位などを持つ。 これを NetCDF というフォーマットを離れ純化したものが GPhys クラスであ る。GPhys データの物理的な所在は、NetCDF などのファイル中にある場合(こ の場合ファイルハンドラーを内部データとして持つ)、メモリー上の配列 (NArray / NArrayMiss) に格納されている場合、そして他の GPhys オブジェ クトのサブセットである場合や、逆に複数の GPhys オブジェクトを束ねる場 合と、多様である。サポートされるファイル形式も現在 3種類ある。これらが 全て、格子点上の物理量として同じ仕様で扱えるようになっている。そもそも、 GPhys クラス自体はこれらの違いを知らず、多様性のサポートはより下位のク ラスに任されている。

物理量は単位を持つため、GPhys は演算時に単位も自動更新する。このため Ruby で書かれた下記の 単位ライブラリー Units を利用する。

サンプル

NetCDF の項で示したサンプルプログラム netcdf_create.rb により作成した ファイルを読み込み、等高線図を表示し、さらに、y 座標に関し -15 度 から +15 度の範囲で平均したものを x の関数として折れ線表示する。ソー スは下のようになる。

gphys_demo.rb

   1|require "numru/ggraph"      # require "numru/gphys" if graphic is not needed
   2|include NumRu
   3|
   4|gphys = GPhys::IO.open('test.nc','u')   # create a GPhys obj on data in file
   5|
   6|DCL.gropn(1)
   7|DCL.sldiv('y',2,1)                   # 'y'oko narabe, 2 by 1 window division
   8|GGraph.contour( gphys )
   9|GGraph.line( gphys.cut('y'=>-15..15).mean('y') )
  10|DCL.grcls

これまで示した可視化プログラムに比べて圧倒的に短い。これを実行すると下 の図が表示される。短いながらも、座標軸まで正しく書いているのがポイント である。

1 行目で読み込んでいる "numru/ggraph" は、GPhys を利用する可視化モジュー ルである(可視化が要らない場合は、GPhys の本体 "numru/gphys" を読み込 めば良い)。2 行目は、これまで同様のおまじないである。そして、4 行目で GPhys オブジェクトを構成している。座標変数は、主変数である 'u' から自 動的に検索されるので指定しなくて良い。6 行目は RubyDCL のグラフィックデ バイス初期化であり、7 行目は画面を 2×1 に分割することを指示する。そし て、8 行目で等高線図を書き、9 行目では "y" という名前の座標に関し、-15 から 15 までの範囲を取出して平均し、折れ線として表示する。9 行目は DCL の 終了処理である。このように、GGraph は DCL の命令と協調して使われるので、 DCL の知識は多少必要である。なお、変数 gphys に割り当てられた GPhys オブ ジェクトは、ファイル中のデータを直接代表しており、数値データは持たない。 GPhys は、巨大なデータに対処できるよう、データ読み込みを実際に必要にな るまで行わない設計になっている。サブセット切り出しも参照のみで行うので、 9 行目で

   gphys.cut('y'=>-15..15)

とした時点でもまだデータはファイル中にあるままである。その次の、平均操 作をする際に初めて、必要な範囲のデータが読み出されるのである。

gphys_demo_sml.png
gphys_demo.rb 実行結果のスクリーンショット

単位ライブラリー Units

Ruby のみで書かれた物理量の単位のクラスである。基本単位への分解機能を 持ち、単位同志の演算や数値の単位変換ができる。 udunits という C のライブラリーに似るが、udunits と違い、演算後も自然な文字列 表現が出来るよう工夫されている。

終わりに

本稿では、多次元数値配列 NArray の紹介を皮切りに、可視化、数値計算、入 出力のライブラリーを紹介した。本人が使っているものを中心とした説明では あったが、Ruby でも大規模な数値データを扱えることは、理解頂けたであろ うか。実際、筆者はスパコンで大気のシミュレーションを行って作った大規模 データの処理に Ruby を使っている(シミュレーションコード自体は Fortran90 で書かれている)。最後に紹介した応用ライブラリーから分かるよ うに、Ruby でライブラリーを作り込めば、数値計算/可視化のプログラミン グの大幅な省力化がはかれる。プログラムが短くなると同時に直感的に分かり やすくなるという恩恵は、少々の遅さを補って余りある効用がある。本当に実 行速度が問題になるところは別の言語で、そうでないところは Ruby でという使 い分けはいかがであろうか。今後もさらに、Ruby の科学計算ライブラリーや 応用例が増えていくことを期待している。

謝辞

本稿を書くよう声をかけてくださった立石孝彰さん、紹介したライブラリーの 作者である田中昌宏さん、西澤誠也さん、常定芳基さん、そして様々な形で電 脳 Ruby プロジェクトに関わられている皆さんに感謝します。

参考資料

その他、本文を参照のこと。

著者について

専門は大気物理学で、京都大学に勤務しています。 学部生の数値解析演習も担当してますが、言語は C です (演習用には、本稿で紹介したような既製のライブラリー は使いません (^_^))。

更新日時:2005/10/02 22:57:19
キーワード:
参照:[Rubyist Magazine 0006 号] [各号目次] [pre-RLRPlan] [prep-0006]