日本語練習帳

つらみのラプソディを歌います

複数のファイル間でデータ平均を取るtips

数値計算による研究をしていると,

  1. 何らかの方法でデータを拾い,無加工で出力する
  2. 出力されたデータを処方箋に則って加工する
  3. 加工されたデータを視覚化し,考察する

というサイクルに耽溺することが多いかと思います.そればかりやるのが良いか悪いかは置いておいても.

そして,多くの場合は出力加工は自動化することが出来て,何も考えなくてもプロセスが進むため,研究している側は他のことに時間を費やすことが出来ます.

しかし,比率で言って出力よりも加工の方が圧倒的に時間が掛かるなと思ったら,何かを疑ったほうが良いでしょう.もしそのような事態に陥ったら,原因として

という可能性を探ってみて下さい.このような組み合わせは,端的に言ってバランスが悪いです.

そこで,思い切ってなるべくスクリプト言語を使わないという方法を提案し,データセットが十分に大きい場合にはこの考えが役立つことを示していこうと思います.


※本記事では実データの例時間計測の結果を示すには時間の都合上至りませんでした.また後日にでも続きを書こうと思います.


始めに

以下のような2つのテキストファイルがあるとします.

1 1.0 1.0
2 1.2 0.9
3 1.4 0.89
4 1.6 0.899
5 1.8 0.8999
1 2.0 2.0
2 2.2 1.9
3 2.4 1.89
4 2.6 1.899
5 2.8 1.8999

そして,それぞれをdata1.datdata2.datと名付け,併せて一つのデータセットだと考えます.

つまり,data{1,2}.datの1-3列目をそれぞれ $$ \begin{align*} t,\quad \{f_{i}(t)\}_{i=1,2},\quad \{g_{i}(t)\}_{i=1,2} \end{align*} $$ とおいて,$f(t)$,$g(t)$という2種類の時間依存データが2つ存在すると見做します1

さて,これらの2つのデータの2列目と3列目の平均を取って,

1 1.5 1.5
2 1.7 1.4
3 1.9 1.49
4 2.1 1.499
5 2.3 1.4999

という新しい1つのデータに統合することを考えます.即ち,数式で表現すると $$ \begin{align*} \tilde{f}^{(N)}(t):=&\frac{1}{N}\sum_{i=1}^{N}f_{i}(t)\\ \tilde{g}^{(N)}(t):=&\frac{1}{N}\sum_{i=1}^{N}g_{i}(t) \end{align*} $$ という新しい量を定義し,それらのデータを $$ \begin{align*} t,\quad \tilde{f}^{(2)}(t),\quad \tilde{g}^{(2)}(t) \end{align*} $$ という形で吐き出すことを目標に据えます.

一般にこういう作業を「クローンデータの平均化」と読んだりします2

こうやって書くと比較的シンプルな課題なので,シェルスクリプトとか単体のawkスクリプトによってさっくりと解決できそうなものです.

実際,データセットのサイズが十分に小規模のものについては,幾つかのブログで紹介されています.自分は,特に『複数ファイルの列データを1ファイルにまとめる』を参考にしたので,挙げておきます.

masa-cbl.hatenadiary.jp

短く言えば,「スクリプト言語の範囲で簡単に出来ますよ」という内容です.


スクリプト言語からの逸脱

ところが,上に挙げたようなデータセットのサイズを遥かに超えるような巨大なものについては,何が起きるかわからないのがスクリプト言語というものです.つまり処理速度の点では「向かない」と言えます.例えば,一般に $$ \begin{align*} t,\quad \tilde{f}^{(N)}(t),\quad \tilde{g}^{(N)}(t) \end{align*} $$ と書いて,$N$が$1000$とか$10000$のオーダーで,$t$方向にも$1000$行~$10000$行程度のデータが並ぶような場合を考えます.そのようなサイズのデータに対して,同様の処理をスクリプト言語でやろうとすると,(データ構造にも依りますが)よっぽど気を使って最適化しない限り1日や2日の規模の時間を消費します3.そして,普通は1日をかけるような贅沢は言っていられないわけです.

やがては「いっその事スクリプト言語を辞めてC++Fortranなどのコンパイル式の言語に頼った方が良いんじゃないの?」という問に到達します.「パンがないならお菓子を食べればいいじゃない」ではないですが,時間との戦いが厳しくなってきたと思ったら,発想の転換を余儀なくされます4


Fortranによるデータセットの平均化

そこで,データ加工をFortran(=Fortran95)を使って行う方法を紹介します.

まず,ソースコードを以下に示します.以下のコードは

  • 加工したいデータセットの一つ一つがdata001.datdata002.dat,…という風にゼロパディングされた3桁の自然数によって番号付けられていて,実行ディレクトリ直下に存在すること
  • それらの個数(n)及び各データの行数(time)をユーザーが把握していること
  • 実行ディレクトリ直下にdata.datというファイルが存在しないこと

を前提として動きます5

PROGRAM analyze

  IMPLICIT NONE

  !変数宣言
  CHARACTER si*4
  INTEGER(kind = 4) :: i, it, dum_t, dum_d
  INTEGER(kind = 4) :: n, time
  DOUBLE PRECISION, ALLOCATABLE :: d(:)

  !標準エラー出力(装置番号0)でコンソール部分を分離する
  WRITE(0, *) "[Input] Number of Clones."
  READ(*, *) n !data*.datの個数を標準入力に読ませる
  WRITE(0, *) "[Input] Length of each Stream."
  READ(*, *) time !各data*.datの行数を標準入力に読ませる

  ALLOCATE(d(1:time)) !データを読むための配列を割付

  d(1:time) = 0d0 !配列を初期化

  !データを配列にスタック
  DO i = 1, n, 1
     WRITE(si, '(I0.3)') i !数値を文字列に変換
     OPEN(10, file="data"//si//".dat", status="old")
     DO it = 1, time, 1
        READ(10, *) dum_t, dum_d !ダミー変数に読み込ませる
        d(it) = d(it) + dum_d !配列に順次加算
     END DO
     CLOSE(10)

     OPEN(11, file="data.dat", status="new")
     DO it = 1, time, 1
        WRITE(11, *) it, d(it)/n !平均化されたデータを出力
     END DO
     CLOSE(11)
  END DO

  DEALLOCATE(d) !配列の割付を解除

END PROGRAM analyze

次に,プログラム内部の流れはコメントから明らかだと思うので割愛することにし,これを使ったり,改変して別の用途に用いたりするための普遍的な手順を紹介します.

  1. 処理したいデータの全てを実行ディレクトリ直下に配置する.
  2. 上記のソースをコンパイルして実行する.
  3. 標準入力からデータの数各データの行数を入力する.

多くの場合,そうやって出てきた加工済みデータ(data.dat)に対して更なる処理を加えたり,gnuplotで視覚化したりすることが想定されるため,加工処理が済んだら,加工前データ(data???.dat)を適当に退避させると良いでしょう.

例えばGNUコンパイラの場合だと

gfortran main.f90 -o main

などで実行ファイルを作ってから

echo "${numParam} ${numClone}" | ./main
mkdir clone
mv data???.dat clone/

を1セットで行うのもありです.…と思った矢先に

/bin/mv: Argument list too long

と怒られる事案が発生しました.こういう時は,まずmvの代わりにxargsを用いて

find . -name "data???.dat" | mv {} clone/

を試すのがお約束でしょう.しかし,それでも怒られることがあります.その時は,最終手段としてGNU Parallelを用いて全CPUを叩き起こせば何とかなると思います.

find . -name "data???.dat" | parallel --bar -a - mv {} clone/

何とかならなかった人は,ご愁傷様です.この優先順位でparallelに手を出した場合、それなりに時間がかかることが予想されるため、一応--barによって進捗を観察できるようにしています.

ひとまず,試してみたい方はお手持ちのコンパイラ6で上記のソースをコンパイルして,爆速のデータ処理を存分にご賞味あれ.


シェルスクリプトによるデータセットの平均化

比較対象と言ってはなんですが,自分で試行錯誤の末に作ったシェルスクリプトワンライナー版を以下に示します7

seq 1 1 $time | parallel -k --bar -a - "awk '\$1=={1} {d+=\$2} END{print {1}, d/NR}' data*.dat" >> data.dat

各々のシェル変数は,Fortranソースコードとの対応を考えて,同じ文字を使っています.nに対応するシェル変数は正規表現のお陰でどこかに飛んでいってくれました.

--barは,何となく「進捗が目に見えないと気が狂いそうだから」という理由で付けています。

これを使ってデータ処理を行うと,飽くまで記憶の範囲ですが,かなりの時間を取られました.「parallelの効果がCPUのコア数に依存するから」というのも考えられる一つの原因ですが,メモリアクセスの効率も一般にC++Fortranに比べて悪い可能性があります.

逆に言えば,そこを改善できれば,処理速度の問題は解決できると踏んでいます.


終わりに

本当は,例えばawkワンライナーC++とかFortranに匹敵するスピードでデータ処理を行うのが目標でした.

しかし,速度比にして10倍~100倍の悦びを知ってしまった今,思いつきの「(大きなデータを処理する時は)スクリプト言語を使うのをやめよう」という考えにシフトしつつあります.

ただ,スクリプト言語についてはとにかく勉強が浅いので,ちゃんとやれば任意のスクリプト言語で爆速データ処理が出来るという希望も捨てずに,色々な可能性を探っていこうと思います.

一方で,この記事の内容は「コマンドの起動時間とデータの参照時間の2つは飽くまで兼ね合いの関係にあるのだから,それはそう」と言われてしまえばお終いです.でも,それを自分の手で確かめることは大いに意義があると思われます.


(2017/02/16追記)

Mathjaxの調子がおかしくて何度も更新を繰り返してしまった.

近いうちに,はてブMarkdown記法で数式を記述する時に必要な手順についてまとめようと思います.


  1. ここでは,両者ともに1列目($1,2,3,4,5$)が共通しているのがポイントで,個々で言う「データ」には情報としては含まれません.逆に,これらが異なっている場合には,シェルスクリプトseqpasteを使って$1,2,\ldots$という縦のストリームを左側に付け加えれば,この記事の議論がそのまま適用できます.

  2. 勘の良い方はお気づきかも知れませんが,自分は物理のシミュレーションをする際に必要に迫られ,このテクニックを取り入れました.例えば,Ising模型における磁化やエネルギーなど,(一般に)MC時間に依存するような物理量の時間発展は,初期状態と乱数の種子によって決まります.つまり,初期状態を決めてシミュレーションを行っても,乱数の種子を入れ替える限り,無数の種類の(系の)時間発展が考えられます.一方,統計力学では,しばしば物理量を「圧倒的に多くの互いに無相関なサンプルの平均」によって決定します.この無相関なデータセットとして「ある初期状態を決めた時の,幾つかの乱数の種子($i=1,2,\ldots$)に対するある量$\hat{F}$の時間発展の組$\{F_{i}(t)\}_{i=1,2,\ldots}$」というものが妥当であれば,それを採用しない手はありません.そこで,これらのデータセットのそれぞれの種子に対するものを「クローン」と呼び,それらを全部平均したものを「クローン平均」と呼ぶことにします.名前の由来は忘れましたが,互いに異なるクローンは(量子力学的な場合を除いて)MCMCの範囲では同時に実現しないからこそ「クローン」なんだろうなと推測しています.もし詳しい方がいらっしゃったら,ご教示願いたいです.

  3. 勿論,こんな主張は要出典なわけで「データを参照する順序がおかしいから時間を消費するんだ」とか「parallelコマンドの使い方を間違えてるんじゃないの」とか,ツッコミの可能性はいくらでもあります.ただ,ここで取り扱いたいのはあまり何も考えずに爆速でデータ処理をするといういわば脳筋的なテーマなので,悪しからず.

  4. ポール・ジョンソンの『インテレクチュアルズ』にもあるように,この言葉にまつわるエピソードの言い出しっぺはルソーで,一方で彼自身は全く出典を示していないそうです.少なくとも,マリー・アントワネットの言葉ではないらみたいで,自分もこの前Twitterで知ったばかりです.

  5. これらの前提条件を崩したい場合には,牛島さんの『数値計算のためのFortran90/95プログラミング入門』などを参照しつつ上のコードを書き換えて下さい.

  6. 大概はGNUコンパイラのgfortranかIntelのParallel Studio XEに付属しているやつのいずれかでしょう.

  7. 無論,ちゃんとやればFortran版に匹敵するスクリプトを組める可能性はある.また,PerlとかRubyとかPythonについては不勉強故コードの書き方を知らなかったため,今回は試していない.案外,そういう所に落とし穴があるのかも知れない.何かアイデアがあれば,ご教示願いたいです.