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

初稿：2005-05-09

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

はじめに

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

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

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

登録データ

RAAリンク: RAA:narray
作者: 田中昌宏
ホームページ: http://www.ir.isas.ac.jp/~masa/ruby/
バージョン: 0.5.7
レポート環境: ruby 1.8, 1.6
必要なライブラリー: (なし)

概要

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 添字が列で、第２添字が行と、行列ならぬ「列行」になっている。よって、本来の行列の意味に立ち返ると、実は列に関して連続で行に関してより「大きく」回る、 row (行)-major であることは注意を要する。このため、添字順に「行列」にする慣例に従っている外部ライブラリーに対しては、むしろ C との相性がよい。

サンプル

初歩

最初に簡単な例を示す。

na1.rb

require "narray"
na = NArray.float(3,2).indgen!         # indgen! --> 0,1,2,..
nb = NArray.to_na( [ [1.0],            # Will be a float NArray
                     [10 ] ] )
nc = NArray.sfloat(12).random!
print "na = "; p na
print "nb = "; p nb
print "nc = "; p nc
print "na + nb = "; p na + nb
print "na >= nb = "; p na >= nb
na[na >= nb] = 999
print "na[na>=nb]=999 --> na = "; p na

print "na[0..1,true] = "; p na[0..1,true]
print "na[true,-1]";      p na[true,-1]
print "na[[0,2,0],true]"; p na[[0,2,0],true]

include NMath
print "sqrt(na) = "; p sqrt(na)
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 行目では、３つの 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 行目では大小比較を行っている。２項演算は、同じ形の配列に対する要素毎の演算となる。しかし、長さ１の次元は任意の長さの次元に拡張されることになっているため、上のように 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

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

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

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

na3.rb

require "narray"

class NArray
  # differentiate along the dim-th dimension:
  #   self[.., 1..-1, ..] - self[.., 0..-2, ..]
  def diff(dim)
    idx1 = [true]*dim + [1..-1, false]
    idx2 = [true]*dim + [0..-2, false]
    self[*idx1].sbt!( self[*idx2] )   # sbt! : destructive subtraction
  end
end

nx = 5
ny = 4
x = NArray.float(nx,1).indgen!
y = NArray.float(1,ny).indgen!
na = x * y
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

require "narray"

def hilbert(n)
  h = NMatrix.float(n,n)
  for i in 0...n
    for j in 0...n
      h[j,i] = 1.0 / (i+j+1)
    end
  end
  h
end

[3,6,12,13].each do |n|
  begin
    h = hilbert(n)
    e = NMatrix.float(n,n).unit
    hi = e / h
    print "n=#{n}. Error in h*h^-1 : ", (e - h*hi).abs.max, "\n"
  rescue
    print "n=#{n}. COMPUTATION FAILED: ",$!,"\n"
  end
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

登録データ

RAAリンク: RAA:rubydcl
作者: 堀之内武ほか
ホームページ:
- RubyDCL: http://ruby.gfd-dennou.org/products/ruby-dcl/
- DCL本家: http://www.gfd-dennou.org/arch/dcl/
バージョン: 1.5.0
レポート環境: ruby 1.8, 1.6
必要なライブラリー: NArray, DCL (Cバージョン)、NArrayMiss

概要

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

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

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

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

サンプル

1次元描画

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

dcl_hop.rb

require "narray"
require "numru/dcl"
include NumRu                 # Most dennou-ruby products are in it
include NMath                 # Math module for NArray

n = 400
dt = 2*PI/(n-1)
t = NArray.sfloat(n).indgen! * dt        # 0 <= t <= 2*pi
x = 1e2 *sin(4*t)
y = 1e-3*cos(5*t)+6

DCL::gropn(1)                      # 1=>Display; 2=>PS file; 4=>Gtk
DCL::grfrm                               # initilize a frame (page)
DCL::ussttl('X-TITLE', 'x-unit', 'Y-TITLE', 'y-unit')
DCL::usgrph(x, y)
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.rb 実行結果のスクリーンショット

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

2次元描画

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

dcl_contour.rb

=begin
=dcl_contour.rb: 等高線と色塗りのサンプル
使用法
 % ruby dcl_contour.rb [-ps] [-dump]
履歴
 * 水田亮 作成 (「ごくらくDCL」デモプログラムを改変)
 * 堀之内 改造 2005/04
=end

require "numru/dcl"
include NumRu
include NMath

nt = 50
nz = 50
tmin, tmax = 0.0, 5.0
zmin, zmax = 20.0,50.0

t = NArray.sfloat(nt+1,    1).indgen! * (tmax-tmin)/nt    # 横軸
z = NArray.sfloat(   1, nz+1).indgen! * (zmax-zmin)/nz    # 縦軸

uz = exp(-0.2*z)*sqrt(z)
tz = -2.0*exp(-0.1*z)
u = uz*sin(3.0*(tz+t))               # (nt+1,nz+1)の配列になる。

if ARGV.index("-ps")
  # コマンドラインオプションに-psが含まれていればPSファイルを作る
  DCL::gropn(2)     # Windows では 3 にすること
else
  # ウィンドウを開く
  DCL::swpset('iheight',700)           # 画面の縦方向のピクセル数指定
  DCL::swpset('iwidth',700)            # 画面の横方向のピクセル数指定
  DCL::gropn(1)     
  DCL::swpset('ldump', true) if ARGV.index("-dump")   # 画面のダンプ
end

DCL::grfrm
DCL::grswnd(tmin, tmax, zmin, zmax)         # 左・右・下・上の座標値
DCL::uspfit                                 # あとはおまかせ
DCL::grstrf                                 # 座標の確定

DCL::ussttl('TIME', 'YEAR', 'HEIGHT', 'km') # 軸のタイトルと単位の設定
DCL::uelset("ltone",true)                   # 色つきにする
DCL::uetone(u)                              # ぬりわけ
DCL::usdaxs                                 # 座標軸を描く
DCL::udcntr(u)                              # 等値線をひく
DCL::uxmttl('t','Idealized QBO',0.0)        # title. top に中央合せ(0.0)

DCL::grcls

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

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

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

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

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


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

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

Ruby-VTK

RAAリンク: RAA:ruby-vtk
作者: 西澤誠也
ホームページ: http://ruby.gfd-dennou.org/products/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 作者が書いたサンプルプログラムを使って、ある数値データの等値面を描画した際のスクリーンショットである。

Ruby/PGPLOT

RAAリンク: RAA:ruby-pgplot
作者: 田中昌宏
ホームページ: http://www.ir.isas.ac.jp/~masa/ruby/#pgplot

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

その他

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

レポート Ruby/GSL

登録データ

RAAリンク: RAA:rb-gsl
作者: 常定芳基
ホームページ:
- Ruby/GSL: http://rb-gsl.rubyforge.org/
- GSL本家: http://www.gnu.org/software/gsl/
バージョン: 1.6.2
レポート環境: ruby 1.8
必要なライブラリー: GSL 及びオプションで NArray

概要

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

require "gsl"

def hilbert(n)
  # 実は GSL::Matrix にはクラスメソッドとして hilbert があるのだがデモのため
  h = GSL::Matrix.new(n,n)
  for i in 0...n
    for j in 0...n
      h[i,j] = 1.0 / (i+j+1)
    end
  end
  h
end

[3,6,12,13].each do |n|
  h = hilbert(n)
  hi = GSL::Linalg::LU.invert(h)
  e = GSL::Matrix.unit(n)
  diff = e - h*hi
  err = [ diff.max, -diff.min ].max      # diff.abs.max if rbgsl >= 1.6.3
  print "n=#{n}. Error in h*h^-1 : ", err, "\n"
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

require "narray"
require "gsl"

def hilbert(n)
  h = NMatrix.float(n,n)
  for i in 0...n
    for j in 0...n
      h[j,i] = 1.0 / (i+j+1)
    end
  end
  h
end

[3,6,12,13].each do |n|
  h = hilbert(n)
  ghi = GSL::Linalg::LU.invert(h.to_gm_view)
  hi = NMatrix.refer( ghi.to_na )       # Simply ghi.to_nm if RubyGSL >= 1.6.3
  e = NMatrix.float(n,n).unit
  print "n=#{n}. Error in h*h^-1 : ", (e - h*hi).abs.max, "\n"
end

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

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

Ruby-SSL2

RAAリンク: RAA:ruby-ssl2
作者: 西澤誠也
ホームページ: http://ruby.gfd-dennou.org/products/ruby-ssl2/

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

Linalg

RAAリンク: RAA:linalg
作者: Jeff Mitchell
ホームページ: http://linalg.rubyforge.org/

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

Linalg の LAPACK モジュールは、FORTRAN77 による LAPACK のオリジナルライブラリーを１対１にラップする。その戦略は、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

登録データ

RAAリンク: RAA:rubynetcdf
作者: 堀之内武ほか
ホームページ:
- RubyNetCDF: http://ruby.gfd-dennou.org/products/ruby-netcdf/
- NetCDF本家: http://www.unidata.ucar.edu/packages/netcdf/, http://my.unidata.ucar.edu/content/software/netcdf/
バージョン: 0.6.0
レポート環境: ruby 1.8, 1.6
必要なライブラリー: NArray, NArrayMiss

概要

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

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

NetCDF のデータモデル

NetCDF の構成を示す。[ ] 内は RubyNetCDF におけるクラス名である。データセット中には、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

require "numru/netcdf"
include NumRu
include NMath

#   -- create a new dataset --
file = NetCDF.create("test.nc")             # .nc is the standard suffix

#   -- define data --
nx, ny = 50, 25
xdim = file.def_dim("x",nx)
ydim = file.def_dim("y",ny)
require "date"
file.put_att("history","Created by #{$0}  #{Date.today}")

xvar = file.def_var("x","sfloat",[xdim])        # x coordinate variable
xvar.put_att("long_name","longitude")
xvar.put_att("units","degrees")

yvar = file.def_var("y","sfloat",[ydim])        # y coordinate variable
yvar.put_att("long_name","latitude")
yvar.put_att("units","degrees")

uvar = file.def_var("u","sfloat",[xdim,ydim])   # 2D along x and y 
uvar.put_att("long_name","wind speed")
uvar.put_att("units","m/s")

file.enddef                                  # end the define mode

#   -- create data values for test --
x = NArray.sfloat(nx, 1).indgen! * (360.0/(nx-1))
y = NArray.sfloat(1, ny).indgen! * (180.0/(ny-1)) - 90.0
x_rad = x * (PI/180)
y_rad = y * (PI/180)
u = (30 + 5*cos(2*x_rad)) * ( 0.2 - cos(3*y_rad) )

#   -- write data values --
xvar.put( x )
yvar.put( y )
uvar.put( u )

file.close

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

イテレーター

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

netcdf_copy.rb

#  Makes a copy of a NetCDF file:
#  % ruby netcdf_copy.rb path_from path_to

require "numru/netcdf"
include NumRu
usage = "\n\nUSAGE:\n% ruby #{$0} path_from path_to\n"
raise usage if ARGV.length != 2
path_from, path_to = ARGV
from = NetCDF.open(path_from)
to = NetCDF.create(path_to)
from.each_dim{|dim| to.def_dim(dim.name, dim.length_ul0)} # len==0 if unlimited
from.each_att{|att| to.put_att(att.name, att.get)}
from.each_var do |var|
  newvar = to.def_var(var.name, var.ntype, var.dim_names)
  var.each_att{|att| newvar.put_att(att.name, att.get)}
end
to.enddef
from.each_var{|var| to.var(var.name).put(var.get)}
to.close

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

読み込み

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

netcdf_plot.rb

require "numru/netcdf"
require "numru/dcl"
include NumRu

#  -- read data --
file = NetCDF.open("test.nc")          # made with netcdf_create.rb
xvar = file.var("x")                   # get a NetCDFVar object named "x"
yvar = file.var("y")
uvar = file.var("u")
x = xvar.get                           # read values of the variable
y = yvar.get
u = uvar.get
xname = xvar.att('long_name').get      # get a NetCDFAtt and read values
yname = yvar.att('long_name').get
uname = uvar.att('long_name').get
xunits = xvar.att('units').get
yunits = yvar.att('units').get
uunits = uvar.att('units').get

# -- graphics --
DCL.gropn(1)
DCL.sgpset('lfull',true)       # use the whole window even if not a square
DCL.uzfact(0.8)                             # --> a bit smaller characters
DCL.grfrm
DCL.grswnd(x.min, x.max, y.min, y.max)        # physical boundaries
DCL.grsvpt(0.2, 0.8, 0.2, 0.6)                # viewport on window
DCL.grstrf                                    # fix the "transform"
DCL.uelset("ltone",true)                      # color shading from now on
DCL.uetone(u)                                 # shading
DCL.ussttl(xname, xunits, yname, yunits)      # set axis titles
DCL.usdaxs                                    # axis
DCL.udcntz(u)                                 # contour
DCL.uxmttl('t', "#{uname} [#{uunits}]", 0.0)  # title at the top
DCL.grcls

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

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

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

レポート GPhys

登録データ

RAAリンク: RAA:gphys
作者: 堀之内武ほか
ホームページ: http://ruby.gfd-dennou.org/products/gphys/
バージョン: 0.4.0
レポート環境: ruby 1.8, 1.6
必要なライブラリー: NArray, NArrayMiss, RubyNetCDF, Units, NumRu::Misc, RubyDCL. オプションで Ruby-FFTW3, Ruby-SSL2.

概要

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

require "numru/ggraph"      # require "numru/gphys" if graphic is not needed
include NumRu

gphys = GPhys::IO.open('test.nc','u')   # create a GPhys obj on data in file

DCL.gropn(1)
DCL.sldiv('y',2,1)                   # 'y'oko narabe, 2 by 1 window division
GGraph.contour( gphys )
GGraph.line( gphys.cut('y'=>-15..15).mean('y') )
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.rb 実行結果のスクリーンショット

単位ライブラリー Units

RAAリンク: http://raa.ruby-lang.org/project/units/
作者: 豊田英司ほか
ホームページ: http://ruby.gfd-dennou.org/products/numru-units/

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

終わりに

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

謝辞

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

参考資料

地球流体電脳 Ruby プロジェクトホームページ: http://ruby.gfd-dennou.org/
- 今回紹介した多くのライブラリーと、rpm, dev や windows 用パッケージなど。
「Ruby for Scientific Computing」: http://pub.cozmixng.org/~the-rwiki/rw-cgi.rb?cmd=view;name=Ruby+for+Science
- Ruby で科学計算するためのライブラリーやプロジェクトへのリンク集。

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

著者について

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

バックナンバー

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

はじめに

始める前に

レポート NArray

登録データ

概要

サンプル

初歩

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

行列計算

レポート RubyDCL

登録データ

概要

サンプル

1次元描画

2次元描画

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

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

Ruby-VTK

Ruby/PGPLOT

その他

レポート Ruby/GSL

登録データ

概要

サンプル

GSL のみ使う場合

GSL と NArray の連携

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

Ruby-SSL2

Linalg

レポート RubyNetCDF

登録データ

概要

NetCDF のデータモデル

サンプル

データ書出し

イテレーター

読み込み

レポート GPhys

登録データ

概要

サンプル

単位ライブラリー Units

終わりに

謝辞

参考資料

著者について

Ruby Library Report 連載一覧

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