pack関数の使い方
Fortranの配列用組み込み関数の中にpack関数というものがあります．直観的に配列を扱うことができ，データをまとめたりする際に大変便利な関数です
pack関数は配列xの要素の内，条件を満たす要素（正確にはmaskで指定した論理型配列の成分が真である部分の要素*1）のみを取り出して新しい配列を作る関数です．このように書いてもわかりにくいので，具体的に使い方を示します．
以下のような魚の調査をしたとしましょう．
- Station 1からStation 10において，ある魚種の体長を調べる
- 各Stationの底質が砂（Sand）か砂利（Gravel）かを記録する
この調査から以下ようなデータ（data.csv）が得られます．なお，全て架空データです．
このデータから，最小体長，最大体長，平均体長，標準偏差を，調査個体全体，砂場の個体のみ，砂利場の個体のみの三つの区分で算出します．
入力データ
Station, River_bed, TL(cm)
1, Sand , 10.0
1, Sand , 17.2
1, Sand , 19.3
1, Sand , 17.4
1, Sand , 18.5
2, Gravel , 13.1
2, Gravel , 18.8
2, Gravel , 18.9
3, Gravel , 11.8
3, Gravel , 19.6
4, Sand , 11.7
4, Sand , 17.3
4, Sand , 14.0
4, Sand , 18.0
5, Sand , 20.1
5, Sand , 22.0
6, Sand , 16.0
7, Gravel , 17.3
7, Gravel , 16.5
7, Gravel , 12.0
8, Sand , 19.8
8, Sand , 16.7
8, Sand , 14.9
9, Gravel , 18.4
9, Gravel , 13.1
10, Sand , 13.6
10, Sand , 11.5
10, Sand , 15.9
プログラム
module stat_calc implicit none contains subroutine input(unit_num, file_name, st, bed, tl) integer, intent(in)::unit_num character(*), intent(in)::file_name integer, allocatable, intent(out)::st(:) character(*), allocatable, intent(out)::bed(:) double precision, allocatable, intent(out)::tl(:) ! Loop integer line_counter integer i open(unit_num, file=file_name, status="old", action="read") line_counter = 0 read(unit_num, '()') ! ヘッダー行読み飛ばし do read(unit_num, *, end=99) line_counter = line_counter + 1 end do 99 continue rewind(unit_num) ! ファイルを開きなおして，一行目から読み直し allocate(st(line_counter)) allocate(bed(line_counter)) allocate(tl(line_counter)) read(unit_num, '()') ! ヘッダー行読み飛ばし do i = 1, line_counter read(unit_num, *) st(i), bed(i), tl(i) end do close(unit_num) end subroutine ! 平均を求める function dmean(x) result(rslt) double precision, intent(in)::x(:) ! Temporary double precision rslt rslt = dble(sum(x)) / dble(size(x)) end function ! 標準偏差を求める function dsd(x) result(rslt) double precision, intent(in)::x(:) ! Temporary double precision var double precision rslt var = dot_product(x - dmean(x), x - dmean(x)) / dble(size(x)) rslt = sqrt(var) end function end module program main use stat_calc implicit none ! Input integer, allocatable::st(:) ! 調査地点番号 character(len=6), allocatable::bed(:) ! 底質 double precision, allocatable::tl(:) ! 全長 ! Temporary double precision, allocatable::sand_tl(:) ! Sandの場のtlだけを集めたもの double precision, allocatable::gravel_tl(:) ! Gravelの場のtlだけを集めたもの integer n_sand ! 砂場の総データ数 integer n_gravel ! 砂利場の総データ数 call input(17, "data.csv", st, bed, tl) ! packした時の大きさを取得 ! pack関数の書き方は，pack(packされるもの，maskの条件式) ! bed=="Sand"により，bedと同じ大きさで，"Sand"の部分に相当する要素が.true.， ! "Gravel"に相当する要素が.false.である論理型配列が指定される ! bed=="Gravel"の場合も同様に論理型配列を指定している n_sand = size(pack(tl, bed=="Sand")) n_gravel = size(pack(tl, bed=="Gravel")) allocate(sand_tl(n_sand)) allocate(gravel_tl(n_gravel)) sand_tl(:) = pack(tl, bed=="Sand") ! Sandの部分のtlを抽出 gravel_tl(:) = pack(tl, bed=="Gravel") ! Gravelの部分のtlを抽出 print*, "砂場のデータ数 ", n_sand print*, "砂利場のデータ数", n_gravel open(17, file="summary.csv", status="replace") write(17, '(3(A, ","), A)') "", "Total", "Sand", "Gravel" write(17, '(A, ",", 2(i3, ","), i3)') "No_individual", size(tl), n_sand, n_gravel write(17, '(A, ",", 2(f5.1, ","), f5.1)') "Min_TL" , minval(tl) , minval(sand_tl), minval(gravel_tl) write(17, '(A, ",", 2(f5.1, ","), f5.1)') "Max_TL" , maxval(tl) , maxval(sand_tl), maxval(gravel_tl) write(17, '(A, ",", 2(f5.2, ","), f5.2)') "Ave_TL" , dmean(tl) , dmean(sand_tl), dmean(gravel_tl) write(17, '(A, ",", 2(f5.2, ","), f5.2)') "SD_TL" , dsd(tl) , dsd(sand_tl), dsd(gravel_tl) close(17) end program
画面出力
砂場のデータ数 18 砂利場のデータ数 10
出力ファイル
, Total, Sand, Gravel
No_individual, 28, 18, 10
Min_TL, 10.0, 10.0, 11.8
Max_TL, 22.0, 22.0, 19.6
Ave_TL, 16.19, 16.33, 15.95
SD_TL, 3.08, 3.14, 2.96
今回のデータでは区分は底質だけですが，実際の調査ではもっと多くの類型にわけられます．例えば，植物の有無や他魚種の有無でデータを類型化し，類型ごとに今回のような基本統計量やヒストグラムの作成をすることがしばしばあります．
このような時にpack関数によっていじりたいデータだけを抽出すると，プログラムが読みやすく，しかも必要に応じて変更しやすいものになります．
*1:maskや論理型についてはこちら y-tis.hatenablog.com
余白を作らずにepsからpngに変換する
計算やデータ解析の結果を図にする時，私はeps形式で出力します．
図を文書に貼る時は，LaTeXならpdfに変換して貼り*1，MS-Wordではpngにして貼ります．
図の変換にはPDF24を使うことが多いのですが，このソフトでepsをpngにするとA4サイズで出力されるため余白ができ，トリミングが面倒です．
特に，複数のグラフを並べる時には，大きさをそろえるのが面倒です．
epsのバウンディングボックスを反映させるか，せめて余白を除いてpngにしたいところです．
そこで，以下のような手続きを踏みます．
これで描画部分以外の余白が取り除かれたpng画像になります．
ただし，この方法はpdfへの変換を経由するため面倒です．他にもよい方法があると思います．
カワヨシノボリという魚
Fortranのことばかり書いているので，今回は生き物について書こうと思います．
今回の主役はカワヨシノボリ（学名Rhinogobius flumineus）です．
ヨシノボリという種類の魚は他にも日本に生息しており，中には一度海に下って再び河川に戻ってくる種類もいますが，カワヨシノボリは一生を河川で過ごします．
河川によってピークは異なると思いますが，5月から7月くらいにかけて繁殖します．以下の写真のような，黄色の卵を石の裏に産み付けます．成熟が進むと目が見えるようになり，卵はすこし褐色を帯びたような色になります．
少し砂に埋まりつつ（はまり石），下に少し隙間が空いているような石をめくると卵が見つかることが多いです．というのも，彼らは石の下を少し掘り起こして隙間を作って産卵します．卵を産んだ後は雄が卵の世話をします．川那部・水野（1989）*1の602，603ページにカワヨシノボリが石の下に穴を掘る過程をとらえたとてもよい写真があるので，興味のある方はぜひ見てみてください．あんな写真を野外で撮ってみたいものです．
仔稚魚は砂州の水際（水深数cm）で見つかります．小さい個体で7 mm位です．分かりにくいですが，以下の写真の真ん中に稚魚がいます．
比較的写真がとりやすい魚だと思います．水中で使えるカメラを初めて使うような方にはおすすめの生き物です．
optional引数と論理型変数を利用して柔軟な関数を作る
配列にデータ（生物の個体数や水位など）を格納して平均を求める場合，doループを使ってもよいですが，Fortranを使うなら配列用の組み込み手続きを活用したいものです．
大きさ10の整数値配列xの成分を
x(:) = (/1, 2, 3, 4, 5, 6, 7, 8, 9, 10/)
で与えます．この時，配列の成分の和を返す関数sumと配列の成分の個数を返す関数sizeを使って，
dble(sum(x)) / dble(size(x))
でxの平均を算出できます．dble()は倍精度実数への変換です．sumは与える配列が実数なら実数型を返しますが，sizeは必ず整数値を返しますので，dbleを忘れずにつけなければなりません．
このような簡単なデータなら特に問題は起こりませんが，
- ある年は個体数調査ができなかった
- 設置していた水位計が流されてしまった
などの理由でデータが欠損している場合もあります．
欠損値に-1をあてることにして，整数値配列yを
y(:) = (/1, 2, 3, 4, 5, -1, 7, 8, 9, -1/)
で与えましょう．-1を避けて平均を算出するとなると，xの場合と同じようには計算できません．
sumにはmaskを指定できます．-1以外の成分の個数を数えるにはcount関数（条件を満たす成分の個数を返す）が使えますので，
dble(sum(y, mask=y > -1)) / dble(count(y > -1)) ! mask= y /= -1, count(y /= -1)でも可
としてyの平均を与えることができます．
以上の計算をまとめて，欠損値がある場合とない場合で平均の計算法を変える関数を定義しましょう．このような場合，optional引数を定義すると使いやすい関数となります．
optional属性を持たせた引数は省略可能になります．欠損値をoptional引数とすることで，欠損値がある時とない時で平均値の演算を使い分ける柔軟な関数になります．書き方は以下通りです．
function iaverage1(x, undef) result(rslt) integer, intent(in)::x(:) integer, optional, intent(in)::undef ! Temporary double precision rslt if (present(undef)) then ! optional引数が存在する時 rslt = dble(sum(x, mask=x > undef)) / dble(count(x > undef)) else ! optional引数が存在しない時 rslt = dble(sum(x)) / dble(size(x)) end if end function
以下注釈です．
- 配列を引数にするときは，x(:)のように大きさをコロンで指定することで，任意の大きさの配列を指定することができます
- intent(in)属性をつけることで，その引数は入力としてのみ使い，値を書き換えたりする用途には使わないことを宣言したことになります（この例だと，x(:)に値を加えたりするとコンパイラがエラーを出してくれるようになり，意図しない演算によるバグを防ぐことができます）
- 関数の名前にaverage1という風にaverageの前にiをつけているのは，引数x(:)が整数型（integer）であることを明示的にするためです
さて，上記の関数ではy < undef）やy /= undefのような場合に対応できません．
- countで指定する条件はsumのmaskで指定する条件と同じ
- countの条件もsumのmaskも論理型配列で指定できる
この二点を利用して以下のようにiaverage1を拡張できます．
function iaverage2(x, mask) result(rslt) integer, intent(in)::x(:) logical, optional, intent(in)::mask(:) ! Temporary double precision rslt if (present(mask)) then rslt = dble(sum(x, mask=mask)) / dble(count(mask)) else rslt = dble(sum(x)) / dble(size(x)) end if end function
とすることで，状況に応じて柔軟に対応できます．
以上のまとめとして，以下のコードを示します．
module stat_calc implicit none contains function iaverage1(x, undef) result(rslt) integer, intent(in)::x(:) integer, optional, intent(in)::undef ! Temporary double precision rslt if (present(undef)) then rslt = dble(sum(x, mask=x > undef)) / dble(count(x > undef)) else rslt = dble(sum(x)) / dble(size(x)) end if end function function iaverage2(x, mask) result(rslt) integer, intent(in)::x(:) logical, optional, intent(in)::mask(:) ! Temporary double precision rslt if (present(mask)) then rslt = dble(sum(x, mask=mask)) / dble(count(mask)) else rslt = dble(sum(x)) / dble(size(x)) end if end function end module program main use stat_calc implicit none ! Loop integer i ! Parameter integer, parameter::n = 10 integer, parameter::x(n) = (/(i, i=1, n)/) integer, parameter::y(n) = (/(i, i=1, 5), -1, (i, i=7, 9), -1/) logical masking(n) print'(A, 10i4)', "x(:) = ", x(:) print'(A, 10i4)', "y(:) = ", y(:) print*, ! sumは配列の各成分の和を返す（mask指定可能）関数 ! sizeは配列の要素数を返す（mask指定不可）関数 ! countは条件を満たす要素数（.true.の個数）を返す関数 print*, "Average of x = ", dble(sum(x)) / dble(size(x)) print*, "Average of y (y > -1) = ", dble(sum(y, mask=y > -1)) / dble(count(y > -1)) print*, ! 上記の計算をiaverage1として定義 print*, "iaverage1(x) = ", iaverage1(x) print*, "iaverage1(y, -1) = ", iaverage1(y, -1) print*, ! countの条件は論理型変数で指定可能 masking(:) = y > -1 print*, "masking(:) = ", masking(:) print*, "Average of y using masking = ", dble(sum(y, mask=masking)) / dble(count(masking)) print*, ! 上記の計算をiaverage2として定義 print*, "iaverage2(x) = ", iaverage2(x) print*, "iaverage2(y, masking) = ", iaverage2(y, masking) end program
結果
x(:) = 1 2 3 4 5 6 7 8 9 10 y(:) = 1 2 3 4 5 -1 7 8 9 -1 Average of x = 5.5000000000000000 Average of y (y > -1) = 4.8750000000000000 iaverage1(x) = 5.5000000000000000 iaverage1(y, -1) = 4.8750000000000000 masking(:) = T T T T T F T T T F Average of y using masking = 4.8750000000000000 iaverage2(x) = 5.5000000000000000 iaverage2(y, masking) = 4.8750000000000000
（2018.07.13追記）
論理型変数についてはこちらの記事を参照してください．
y-tis.hatenablog.com
論理型（logical）変数の使い方
Fortranでは変数の型として
- 整数型 integer
- 実数型 real
- 倍精度実数型 double precision
- 複素数型 complex
- 倍精度複素数型 complex(kind(0d0))
- 論理型 logical
- 文字型 character
が使えます．今日の話題は論理型（logical）変数です．
長い間これはどんな場合に使う変数なのか理解できませんでした．簡単な計算ばかりしているときにはありがたみが全く分からなかったのです．1と0で代用できるような場面ばかりだとまず使うことがありませんでした．
データの加工などをするときにはpack関数などのmaskを伴う配列関数をよく使います．このmaskの指定に論理型が使えます．統計処理をするようになってからこの手の関数を使うことが多くなり，やっと論理型変数のありがたみが分かりました．
かんたんな例を挙げましょう．細かい説明はコメントアウトを読んでください．
program main implicit none ! Loop integer i ! Parameter integer, parameter::n = 10 integer, parameter::x(n) = (/(i, i=1, n)/) integer, parameter::y(n) = (/1, 200, 3, 400, 5, 600, 7, 800, 9, 1000/) logical masking(n) ! 論理型の値は.true.または.false. ! true, falseの前後のピリオドを忘れないように masking(:) = .true. print*, masking(:) masking(:) = .false. print*, masking(:) print*, ! pack関数はこのように使う ! 材料となる配列xの成分の内，mask=で指定した条件がtrueの部分に ! 相当する成分だけをまとめた配列を作る ! 従って，maskの条件はxに関する条件でなくてもよい（以下の配列yによるmaskがその例） print'(A31, 5(" "), 10i3)', "pack(x, mask=mod(x(:), 2) == 1)", pack(x, mask=mod(x(:), 2) == 1) print'(A31, 5(" "), 10i3)', "pack(x, mask=y(:) < 10)", pack(x, mask=y(:) < 10) ! pack関数のmask=は省略可能 print'(A31, 5(" "), 10i3)', "pack(x, mod(x(:), 2) == 1)", pack(x, mod(x(:), 2)==1) print'(A31, 5(" "), 10i3)', "pack(x, y(:) < 10)", pack(x, mask=y(:) < 10) print*, ! maskやif文などで書く論理式は.true.または.false.を返すので，実はprint文でこのように出力できる ! なお，以下では配列であることを明示的にするためx(:)と書いている ! それぞれ，mod(x, 2) == 1，x < 5，x >= 5と書いてもよい print'(A17, 5(" "), 10L3)', "mod(x(:), 2) == 1", mod(x(:), 2) == 1 print'(A17, 5(" "), 10L3)', "x(:) < 5", x(:) < 5 print'(A17, 5(" "), 10L3)', "x(:) >= 5", x(:) >= 5 print*, ! pack関数はmask=の値が.true.の部分の要素を集めるので，論理型配列で指定可能 masking(:) = mod(x, 2) == 1 print'(A21, 5(" "), 10i3)', "x(:)", x(:) print'(A21, 5(" "), 10L3)', "masking(:)", masking(:) print'(A21, 5(" "), 10i3)', "pack(x, mask=masking)", pack(x, mask=masking) print*, masking(:) = y(:) < 10 print'(A21, 5(" "), 10i3)', "x(:)", x(:) print'(A21, 5(" "), 10L3)', "masking(:)", masking(:) print'(A21, 5(" "), 10i3)', "pack(x, mask=masking)", pack(x, mask=masking) end program
出力結果
T T T T T T T T T T
F F F F F F F F F F
pack(x, mask=mod(x(:), 2) == 1) 1 3 5 7 9
pack(x, mask=y(:) < 10) 1 3 5 7 9
pack(x, mod(x(:), 2) == 1) 1 3 5 7 9
pack(x, y(:) < 10) 1 3 5 7 9
mod(x(:), 2) == 1 T F T F T F T F T F
x(:) < 5 T T T T F F F F F F
x(:) >= 5 F F F F T T T T T T
x(:) 1 2 3 4 5 6 7 8 9 10
masking(:) T F T F T F T F T F
pack(x, mask=masking) 1 3 5 7 9
x(:) 1 2 3 4 5 6 7 8 9 10
masking(:) T F T F T F T F T F
pack(x, mask=masking) 1 3 5 7 9
以上の例のように，trueかfalseをもとにmaskするというのが重要であり，便利なところです．条件によってmaskの内容が異なる場合には，このような論理型変数が威力を発揮します．
（2018.07.13追記）
論理型変数を利用して関数を作った例はこちら．
y-tis.hatenablog.com
（2018.08.11追記）
pack関数を使うときにも同様の考え方をします．
y-tis.hatenablog.com
ifとcaseと分岐処理削除
Fortranでは分岐処理にifとcaseがあります．条件分岐が三つ以上ある場合，if文では条件を満たすまで上から順にelse ifをたどるのに対し，case文では値の評価をした後，条件を満たすところにジャンプするため，分岐が多い場合にはcaseの方が速いそうです．私は統計処理をする際，文字列を用いた分岐に対してはcaseを使うことが多いです（一文字目が"A"の時に処理aをするなど）．
チューニングの際は，まず条件分岐自体を無くせないかを考えます．意外とifやcaseを使わなくても書けるものです．やたらに条件分岐が多いコードを書いている間は，アルゴリズム自体の理解が浅いことが多いです．
今回はdo loopのカウンターの値に応じて処理を変える場合の計算速度をif文，case文，条件分岐を使わない場合で比較します．実用上このような場面は多く，例えば，
- 移流方程式を時間幅dt=0.1でt=0からt=10まで解き，t=1ごとに解を出力したいとき
- マルコフ連鎖モンテカルロ法で10回ごとにサンプリングを行いたいとき
などです．
以下がソースコードです．
program main implicit none ! Parameter integer, parameter::iter = 10**8 integer, parameter::thin1 = 10 integer, parameter::thin2 = 5 integer, parameter::thin3 = 2 integer, parameter::thin4 = 1 ! Loop integer i, j ! Temporary integer s integer m ! Time integer t0 integer t1 integer t_rate integer t_max s = 0 call system_clock(t0) do i = 1, iter m = mod(i, thin1) if (m == 0) then s = 0 else if (m == thin4) then s = s + thin4 else if (m == thin3) then s = s + thin3 else if (m == thin2) then s = s + thin2 else s = s * 2 end if end do call system_clock(t1, t_rate, t_max) print*, calc_time(t0, t1, t_rate, t_max) s = 0 call system_clock(t0) do i = 1, iter m = mod(i, thin1) select case(m) case (0) s = 0 case (thin4) s = s + thin4 case (thin3) s = s + thin3 case (thin2) s = s + thin2 case default s = s * 2 end select end do call system_clock(t1, t_rate, t_max) print*, calc_time(t0, t1, t_rate, t_max) s = 0 i = 1 call system_clock(t0) do while (i <= iter) ! mod(i, thin1) = thin4 = 1 s = s + thin4 i = i + 1 ! mod(i, thin1) = thin3 = 2 s = s + thin3 i = i + 1 ! mod(i, thin1) = 3--4 do j = thin3+1, thin2-1 s = s * 2 i = i + 1 end do ! mod(i, thin1) = thin2 = 5 s = s + thin2 i = i + 1 ! mod(i, thin1) = 6--9 do j = thin2+1, thin1-1 s = s * 2 i = i + 1 end do ! mod(i, thin1) = 0 s = 0 i = i + 1 end do call system_clock(t1, t_rate, t_max) print*, calc_time(t0, t1, t_rate, t_max) contains function calc_time(t0, t1, t_rate, t_max) result(rslt) integer, intent(in)::t0 integer, intent(in)::t1 integer, intent(in)::t_rate integer, intent(in)::t_max ! Temporary double precision rslt integer diff if (t1 < t0) then diff = t_max - t0 + t1 + 1 else diff = t1 - t0 end if rslt = dble(diff)/dble(t_rate) end function end program
結果
0.39000000000000001 0.34399999999999997 0.18800000000000000
若干caseの方が速いですが，あまり変わりませんね．もう少し差が出ると思っていました．条件分岐がない場合は二倍くらいの速度で計算が終わりました．
計算結果の出力を伴う場合，どれくらい差が出るでしょうか．計算部分を以下のように書き直します．三つ目の計算では，write文を書く位置はiのカウントを増やす前であることに注意しましょう．
open(17, file="ex.csv", status="replace") s = 0 call system_clock(t0) do i = 1, iter m = mod(i, thin1) if (m == 0) then s = 0 write(17, *) i, ",", s else if (m == thin4) then s = s + thin4 write(17, *) i, ",", s else if (m == thin3) then s = s + thin3 write(17, *) i, ",", s else if (m == thin2) then s = s + thin2 write(17, *) i, ",", s else s = s * 2 end if end do call system_clock(t1, t_rate, t_max) print*, calc_time(t0, t1, t_rate, t_max) write(17, *) s = 0 call system_clock(t0) do i = 1, iter m = mod(i, thin1) select case(m) case (0) s = 0 write(17, *) i, ",", s case (thin4) s = s + thin4 write(17, *) i, ",", s case (thin3) s = s + thin3 write(17, *) i, ",", s case (thin2) s = s + thin2 write(17, *) i, ",", s case default s = s * 2 end select end do call system_clock(t1, t_rate, t_max) print*, calc_time(t0, t1, t_rate, t_max) write(17, *) s = 0 i = 1 call system_clock(t0) do while (i <= iter) ! mod(i, thin1) = thin4 = 1 s = s + thin4 write(17, *) i, ",", s i = i + 1 ! mod(i, thin1) = thin3 = 2 s = s + thin3 write(17, *) i, ",", s i = i + 1 ! mod(i, thin1) = 3--4 do j = thin3+1, thin2-1 s = s * 2 i = i + 1 end do ! mod(i, thin1) = thin2 = 5 s = s + thin2 write(17, *) i, ",", s i = i + 1 ! mod(i, thin1) = 6--9 do j = thin2+1, thin1-1 s = s * 2 i = i + 1 end do ! mod(i, thin1) = 0 s = 0 write(17, *) i, ",", s i = i + 1 end do call system_clock(t1, t_rate, t_max) print*, calc_time(t0, t1, t_rate, t_max) close(17)
結果
0.93700000000000006 0.93799999999999994 0.92200000000000004
理由はよくわかりませんが，三者ともほとんど同じ速度になりました．条件分岐のない計算が遅くなったのはなぜでしょうか．今回はここまでとします．
