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
