2024年10月13日日曜日

化石言語FortranをMCMCするのに使ってみたよ!

このエントリーをはてなブックマークに追加
Pocket

数値計算では今でも標準的な地位にあるプログラミング言語Fortranは、最古の高級言語だけに古臭い面が残されており、利用者からの評判も芳しくはない*1。計算用途以外では低機能なので利用者数が伸びるはずもなく、近年、相対的な地位は低下の一方であった。ところがコミュニティの活発さを測るTIOBE Rankingで、2021年頃から順位が上がっている。パッケージマネージャ(fpm)が出てきた影響ではないかと思うが*2、国内外で疑問に思っている人は多いようだ。

Fortranが絶滅予定と言うのは偏見で、私が知らない長所があるのか。長所と短所を体感してみるために、MCMCの一種のHMC法で、線形回帰、混合ロジットモデル、順序ロジットモデル、ポアソン回帰のベイズ統計学による推定と、制約付きモデルのベイズファクターと、ラプラス近似による周辺尤度の計算をするコードを、Fortranで書いてみた。その感想を記しておきたい。

2. 長所

利点は明確にある。

実行すれば高速計算
RのHMCの教育用パッケージhmclearnと比較して、30倍から50倍ぐらい速くなった。科学技術計算をさせればCよりも速いと言われる。これは実際、以前マイクロベンチマークをとったら速かった。
コンパイル時間が短い
Cと同様に速く、C++ほど待たされない。Makefileで分割コンパイルしているからかもだが、Juliaの実行時コンパイルよりも高速かも知れない。
ベクトル演算ができる
RやJuliaに慣れていると当然のように感じるが、CではできないしC++では標準ではない。/(1.0d0, 2.0d0, 3.0d0)/とリテラルの配列をプロシージャに渡すこともできる。行列積をとる関数も標準であり、これも意外に便利である。
並列演算が容易
OpenMPやOpenMPIが使えると言うだけだが、RやPythonでは利用できないので利点となる。GPUサポートもあるが、こちらは試していない。なお、スレッドやセマフォがサポートされているわけではなく、計算以外での並列処理には向かない。
数値計算のためのライブラリは豊富
そもそもFortranで書かれているBLASやLAPACKの他、最適化アルゴリズムもL-BFGS-Bのライブラリがあった。乱数生成はC言語のdSFMTを呼び出すInterfaceが提供されている。
RやPythonからFortranのサブルーチンを呼べる
計算は遅いが便利なインタープリッタに高速なライブラリを提供することができる*3。Cで書くほうが親和性がより高いのでFortranの利点として挙げられない気がするが、Fortranの有力な代替肢となるJuliaはこれが現実的ではない。オブジェクトファイルを生成する設計になっていないからだ。
ちゃんとオブジェクト指向ができる
Fortran 2003からスーパータイプを継承してサブタイプをつくって、スーパータイプを処理するライブラリになげるオブジェクト指向プログラミングができる。今回で言うと、MCMCの部分をフレームワーク化し、それに入れる尤度関数の引数の型をチェックを働かせることができた。

スパコンでのサポートも手厚いと言うような利点もあるそうだ。物理学徒にはうれしいかも知れない。

3. 短所

欠点も明確にある。

ファイル入出力が弱い
ファイル入出力も番号指定で開いたり閉じたり、エラー時はラベル(というか行番号?)に飛ぶ仕様である。I/Oエラーのハンドリングが苦しい。40年前の世界があった。HDF5型式の利用が推奨されているようだ。ただし、データ量が決まっていたらバイナリファイルの読み書きはそう苦にならない。なお、今回はRから呼び出したので、この問題は回避できている。
テキスト処理が苦行
固定長文字列の時代の影響が色濃く、入出力と文字列処理が現在の可変長時代のシステムとあっていない。C言語のprintfに慣れたものは、writeやprintfの書式指定だけで嫌になる。さらに、変数宣言ブロックで固定長の領域を確保して、そこに書式付き文字出力をしてエラーメッセージを作るような煩雑な作業が生じる。
柔軟ではない変数宣言
テキスト処理の問題と重複なのだが、変数宣言や領域の割当に柔軟性がない。宣言文もちょっと冗長である。サブルーチンで配列に領域を割り当てて呼び出し元に返すようなことはできない。構造体(type)が宣言できるのだが、たぶん(C互換型を使わなければ)可変長の要素は持てない。ループ中にループ回数の制御に使っている変数はいじれない。ただし、これらは速度や安全性の都合と言う面もある。
使ってはいけない命令が山ほどある
後方互換性のために1977年規格の文法が使えるが、山ほど問題があるので新たには使ってはいけない命令として非推奨になっている。Modern Fortranと呼ばれる1990年以降の規格で書けと言う話になるのだが、油断していると非Modernになっていることがある。
配列の長さはチェックされない
C/C++と同様ではあるが、コンパイラは配列の添字の範囲が適切かは調べず、プログラマが誤ると謎の挙動を引き起こしたりする。RやPythonやJuliaに慣れていると躓くかも知れない。
型指定が煩雑
整数と単精度浮動小数点と倍精度浮動小数点の自動判別が弱く、リテラルで値を渡すときに明確に指定がいる。0と0e0と0d0の使い分けがいる。
ライブラリの設計が古い
言語仕様上、引数が多くなる傾向があるが、それだけではない。プロシージャの名前が6文字以内というFORTRAN 77時代からのライブラリが現役なので、名前から機能を連想するのが困難だったり、ライブラリが内部で使う作業領域を呼び出し側が確保する仕様となっている。組み込み総称関数が型を自動判別するのに、単精度と倍精度で異なるプロシージャ名を指定するチグハグ感。
NaNやInfが定義されていない
a == aがFALSEになるかでNaNか判定できると言うのをはじめて知った。IEEE 754で規定があると思うのだが。本当に数値計算に特化しているのか。煩雑なので踏襲しなかったが、訓練されたFortran使いは、NaNやInfのIEEE 754定義を参照して、自ら定数定義を与えて使っているようだ。なお、Fortran 2003からieee_arithmeticモジュールが標準になり、そこにはNaNとInfをハンドリングする関数が用意されている。
モジュール・インターフェイス・ファイルがアレ
C言語のインクルードファイルに相当するモジュール・インターフェイス・ファイル(.mod)は自動生成となっている。これが意外に煩雑だ。利用しているモジュールのインターフェイス・ファイルがないとコンパイルができない。ファイル間で中にあるモジュールを相互参照するコードは書くことができない。
ラムダ式が書けない
最近のプログラミング言語に慣れている人は、発想の転換がいる。私はラムダ式や関数ポインタがなくても最適化アルゴリズムのライブラリは作れることを認識した。求めるパラメーターの値をxとして、x₀ → x₁ → … → xₜ → … と漸近するアルゴリズムなのだから、xₜを入れたらxₜ₊₁を戻すプロシージャがあれば足りる。言われてみれば当たり前だが、ループ部分を自分で書くライブラリに前近代を感じた。
方言がある
今回はGFortranだけを使っているので、拡張子.f90と.F90で意味が異なり、後者ではC互換プリプロセッサが標準で有効になることに注意することになったこと以外は顕在化しなかったのだが、以前、Intel FortranのコードをGFortranにポーティングしたときにぽちぽち感じた。C++でもMicrosoft方言があったりするので、Fortranだけの問題ではないが。

今回の作業では大きな問題は無かったが、数値解析以外に使える範囲は極端に狭いと言わざるをえない。多用途(versatile)なにそれ美味しいの? — と言う割り切りは嫌いではないが。なお、スーパータイプを継承してサブタイプをつくって、スーパータイプを処理するライブラリになげるオブジェクト指向プログラミングはできない…と思っていたのだが、このエントリーを書いた後に出来ることに気づいた。

3. 感想

Fortran 90以降(というか2003あたり?)を計算用途で利用したので不自由と言うほどではなかったが、やはり現代的ではない。後方互換性を考えるとこうなってしまう感もあって、新たな科学技術計算用のプログラミング言語が求められているのは分かる。しかし、見よう見まねで計算できるコードが書けてしまっているので、取っ付きにくいかと言うと、そうでもないかも知れない。なかなか新たなプログラミング言語が代替するのも困難そうだ。

*1Fortranではなく他のプログラミング言語を推奨したり、Fortranの利点を示しつつも新たなプロジェクトでの利用はされないと考える人が多い(5 Reasons Why Fortran is Still Used)。少なくともFORTRAN 77は使うな/教えるな、せめてModern Fortranにしろと言われている。

*2チュートリアルの検索回数によるランキングであるPYPL PopularitY of Programming Languageではランク外である。つまり、利用者が増加しているとは考えづらい。

*3RからFortranのサブルーチンを呼ぼう

0 コメント:

コメントを投稿