匿名希望のおでんFortranツヴァイさん太郎

生き物、Fortran、川について書く

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